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A.  SCIENTIFIC  OBJECTIVES:  The  objective  of  this  research  program  is  to  explore 
new  design  methodologies  and  establish  performance  bounds  in  shipboard  antenna 
problems.  Oiu:  approach  is  to  combine  genetic  algorithms  (GA)  and  model-based  signal 
processing  with  computational  electromagnetic  simulators  to  investigate  antenna  design, 
placement  and  optimization.  The  specific  problems  to  be  investigated  include  microstrip 
antenna  design,  array  synthesis  and  placement  in  the  presence  of  platform  effects, 
electrically  small  wire  antenna  design  and  microwave  absorber  design.  Our  goals  are  to 
map  out  fimdamental  perfomumce  bounds  in  complex  design  problems  and  to  devise 
actual  design  implementations  that  can  approach  these  performance  bounds. 

B.  SUMMARY  OF  RESULTS  AND  SIGNIFICANT  ACCOMPLISHMENTS: 
During  the  current  year,  we  have  made  significant  progress  on  four  topics:  (i)  shape 
optimization  of  broadband  and  multi-band  microstrip  antenna  elements,  (ii)  antenna  array 
beamforming  and  placement  in  the  presence  of  platform  effects,  (iii)  shape  optimization 
of  light-weight  microwave  absorbers,  and  (iv)  design  of  electrically  small  wire  antennas. 

In  the  first  topic,  we  have  applied  GA  and  a  fast  computational  electromagnetic  solver 
to  design  novel  microstrip  antenna  shapes  for  broadband  operations,  multi-band 
operations  and  circular  polarization.  For  the  broadband  design,  we  reported  a  four-fold 
inprovement  in  bandwidth  compared  to  a  regular-shaped  microstrip.  For  dual-band 
applications,  we  achieved  arbitrary  fi’equency  spacing  between  the  two  fi’equency  bands 


ranging  from  1.1  to  2.  Tri-band  and  quad-band  operations  were  also  realized.  All 
designs  were  febricated  and  verified  by  measurements. 

In  the  second  topic,  we  have  applied  GA  for  array  beamforming  in  the  presence  of 
platform  and  mutual  coupling  effects.  In  our  optimization  process,  the  active  element 
patterns  of  the  antenna  elements  were  either  computed  using  an  electromagnetic  solver  or 
obtained  through  in-situ  measurements.  The  element  excitations  were  then  optimized  by 
GA  for  each  prescribed  beam.  Both  simulation  results  and  field  measurement  data  (using 
a  7-element  smart  antenna  testbed)  were  collected.  It  was  demonstrated  that  the  GA- 
synthesized  beams  could  be  used  to  overcome  the  blockage  effect  due  the  moimting 
tower.  Direction  finding  experiments  were  also  carried  out  and  showed  that  the  GA 
results  can  be  successfully  used  to  locate  emitters. 

In  the  third  topic,  we  have  applied  GA  to  the  design  of  corrugated  microwave 
absorbers.  By  combining  GA  with  a  fuU-wave  con5)utational  electromagnetic  solver,  we 
obtained  absorber  performance  that  is  superior  to  previously  investigated  geometrical 
shapes.  Physical  interpretation  of  the  GA-optimize  shapes  was  formulated  to  ejqplain  the 
behavior  of  the  absorber  as  a  function  of  polarization  and  incident  angle.  Furthermore, 
we  investigated  the  optimal  absorber  performance  as  a  function  of  absorber  thickness  (or 
weight)  by  means  of  the  Pareto  GA  technique.  The  performance  limit  of  the  absorber  as 
a  function  of  thickness  was  efficiently  mapped  out  using  this  approach. 

In  the  fourth  topic,  we  have  carried  out  research  into  the  design  of  electrically  small 
wire  antennas.  The  well-known  theoretical  limit  for  antenna  bandwidth  as  a  function  of 
its  electrical  size  was  derived  by  L.  J  Chu  in  1948.  Our  objective  was  to  apply  GA 
optimization  to  design  wire  antennas  that  can  approach  the  Chu  limit.  In  om  initial 
design,  we  used  GA  together  with  a  wire  antenna  code  to  study  the  performance  of  a 
smgle-arm  wire  with  a  number  of  bent  segments.  It  was  found  that  the  achievable 
performance  was  below  the  Chu  limit.  Folded  wire  configurations  and  multiple  arm 
designs  are  under  investigation. 

Below  we  shall  describe  the  four  topics  in  more  detail. 

Microstrip  Antenna  Design  and  Optimization.  Due  to  its  low  profile  and  ease  of 
fabrication,  microstrip  antenna  is  a  very  popular  choice  in  many  commercial  and  military 


applications.  However,  a  well-known  drawback  of  the  microstrip  is  that  it  is  intrinsically 
a  narrow-band  device.  We  have  appUed  GA  and  a  fast  computational  electromagnetics 
solver  to  design  novel  microstrip  antenna  sh^es  for  broadband  and  multi-band 
operations  [14,22,24,26].  Research  in  the  application  of  GA  for  antenna  design  has  been 
ongoing  since  the  early  1990s  [1,2].  In  contrast  to  a  local  optimization  algorithm,  GA 
allows  the  design  space  to  be  more  fiilly  explored  (at  the  expense  of  computational  cost). 
In  our  GA  implementation,  the  geometry  of  the  patch  is  discretized  into  a  two- 
dimensional  binary  map  and  GA  is  used  to  search  for  the  optimal  patch  shape  that 
maximizes  bandwidth.  Several  schemes  are  utilized  to  accelerate  the  convergence  of  the 
GA  including  2-point  crossover  and  geometrical  filtering  through  the  use  of  a  median 
filter.  Shapes  in  the  population  are  ranked  according  to  a  computed  cost  function  and  a 
new  generation  of  children  shapes  is  produced  through  the  rules  of  heredity  and  mutation. 
This  process  is  repeated  until  a  satisfactory  shape  is  found  that  best  meets  the  design 
criteria. 

Fig.  1(a)  shows  the  GA-designed  broadband  microstrip  shape  [14].  This  optimized 
shape  has  a  four-fold  improvement  in  bandwidth  when  compared  to  a  standard  square 
microstrip  antenna.  This  design  has  been  verified  by  laboratory  measurement  on  FR-4 
substrate.  Fig.  1(b)  shows  the  excellent  agreement  between  the  simulation  and  measured 
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Fig.  1.  A  broadband  microstrip  antenna  designed  using  the  genetic  algorithm  and  built  on  FR-4 
material.  The  right  figure  shows  the  return  loss  (dB)  of  the  antenna  from  simulation  (□)  and 
measurement  (-).  The  antenna  has  an  operating  fi’equency  of  2  GHz  and  an  8%  bandwidth,  which 
is  4  times  broader  than  a  conventional  square  microstrip. 
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Fig.  2.  Shapes  of  three  GA-optimized  dual-band  microstrip  antennas,  and  the  resulting 

return  loss  from  simulation  ( - )  and  measurement  ( - ).  (a)  Frequency  ratio  of  1:1.3 

(1 .8GHz  and  2.34GHz).  (b)  Frequency  ratio  1:1.6  (1 .8GHz  and  2.9GHz).  (c)  Frequency 
ratio  of  1:1.9  (1.8GHz  and  3.42GHz). 


results.  The  basic  operating  principle  of  the  optimized  shape  has  also  been  interpreted  in 
terms  of  a  combination  of  dual-mode  operation  and  ragged  edge  shape.  We  have  also 
extended  this  methodology  to  design  multi-band  microstrip  antennas  [22,26].  We  have 
designed  a  collection  of  microstrip  shapes  that  can  achieve  two  frequency  bands  of 
operation  with  arbitrary  frequency  spacing  between  the  two  bands  ranging  from  1:1.1  to 
1:2.0.  Figs.  2(a)  to  2(c)  show  three  of  our  designed  shapes  with  frequency  ratios  of  1.3, 
1.6  and  1.9,  respectively.  Again  these  results  have  been  verified  by  measurements  and 
showed  precise  dual-band  operation  with  good  bandwidth.  We  have  also  numerically 
verified  that  these  shapes  could  be  scaled  in  size  to  any  specific  operating  frequency  of 
interest,  or  for  different  substrate  materials.  Tri-band  and  quad-band  microstrips 
designed  using  the  same  GA  methodology  have  also  been  demonstrated  recently  [24]. 

Array  Antenna  Placement  and  Beamforming.  It  is  well  known  in  antenna  design 
that  the  mounting  platform  can  have  a  significant  impact  on  antenna  performance.  For 
example,  in  the  deployment  of  antenna  arrays  it  is  not  always  possible  to  mount  the  array 
on  top  of  a  tower.  When  the  array  is  mounted  at  the  mid-tower  level,  the  pattern  behavior 
of  the  array  could  be  significantly  degraded  due  to  the  presence  of  the  tower.  In  this 
work,  we  address  how  GA  can  be  utilized  to  achieving  array  placement  and  beamforming 
to  overcome  platform  effects.  This  work  is  leveraged  upon  our  previous  studies  on  array 
mutual  coupling  and  platform  interactions  [15,18,19,25]. 

In  our  GA-based  beamforming  procedure,  the  active  element  patterns  of  the  antenna 
elements  are  either  computed  using  an  electromagnetic  solver  or  obtained  through  in-situ 
measurements.  The  element  excitations  are  then  optimized  by  the  GA  for  each 
prescribed  beam.  For  each  beam  pattern,  a  cost  is  assessed  if  the  radiated  field  is  below  a 
desired  level  in  the  main  beam  or  above  the  sidelobe  level  in  the  sidelobe  range.  Unlike 
conventional  array  synthesis  [3]  in  which  the  transition  between  the  main  beam  and 
sidelobe  is  unconstrained,  this  region  is  also  constrained  here  due  to  the  strong  tower 
effect. 

The  algorithm  has  been  tested  using  simulation  data  from  the  Numerical 
Electromagnetics  Code  (NEC)  for  a  seven-element  circular  dipole  array  of  one 
wavelength  diameter.  The  array  is  mounted  next  to  a  metal  tower,  as  shown  in  Fig.  3(a). 


(a)  NEC  simulation 
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Fig.  3.  GA-based  beamforming  by  using  NEC-generated  active  element  patterns. 

(a)  Geometry  of  the  array  with  tower,  (b)  Beam  pattern  of  a  free-standing  array  using 
cophasal  excitation,  (c)  Beam  pattern  of  the  array  next  to  the  tower  using  cophasal 
excitation,  (d)  Beam  pattern  of  the  array  next  to  the  tower  using  GA-optimized  excitation. 
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Fig.  3.  GA-based  beamforming  by  using  measured  active  element  patterns  of  a  7-element 
smart  antenna  array,  (a)  Measurement  setup  of  the  array  with  tower,  (b)  Beam  pattern  of 
the  array  at  150®  using  cophasal  excitation,  (d)  Beam  pattern  of  the  array  using  GA- 
optimized  excitation. 


The  tower  has  a  triangular  cross  section  with  3  wavelengths  on  each  side.  If  the  array  is 
located  in  free  space,  a  beam  can  be  trivially  synthesized  in  any  direction  based  on  the 
cophasal  excitation.  A  sample  beam  is  shown  in  Fig.  3(b).  However,  when  the  tower  is 
present,  the  beam  pattern  degrades  if  the  same  cophasal  excitation  for  the  free-standing 
array  is  used.  This  is  shown  in  Fig.  3(c),  in  which  we  observe  a  degraded  beam  with 
higher  sidelobes.  Finally,  when  GA  is  used  to  synthesize  the  beam,  the  beam  can  be 
nearly  restored  to  the  original  free-space  beam  pattern,  as  shown  in  Fig.  3(d). 

The  algorithm  has  also  been  tested  using  the  measured  active  element  patterns  of  an 
existing  7-element  smart  antenna  array  at  the  University  of  Texas  at  Austin.  The 
measurement  setup  is  shown  in  Fig.  4(a),  where  the  array  is  placed  close  to  a  metal  tower. 
Fig.  4(b)  shows  the  beam  pattern  resulting  from  using  the  standard  cophasal  excitation.  It 
shows  the  effect  of  the  tower  when  the  beam  is  steered  close  to  the  direction  of  the  tower. 
In  particular,  we  see  that  the  beam  is  squinted  away  from  the  desired  150°  direction. 
After  GA  is  performed  to  generate  the  optimal  excitations,  the  resulting  beam  pattern  in 
Fig.  4(c)  is  very  close  to  the  prescribed  beam  pattern.  Thus  the  tower  effect  can  be 
largely  compensated  by  using  the  GA-synthesized  excitation.  Direction  finding 
experiments  have  also  carried  out  and  showed  that  the  GA  results  could  be  successfiilly 
used  to  locate  emitters  in  direction  finding  applications. 

Light-Weight  Microwave  Absorber  Design.  We  have  applied  GA  to  the  design  of 
corrugated  microwave  absorbers  for  oblique  incidence  [27]  (see  Fig.  5).  Previously,  GA 
had  been  applied  to  the  design  of  multi-layered  planar  and  cylindrical  absorber  structures 
[4].  Corrugated  coatings  with  non-planar  shape  profile  offer  an  additional  degree  of 
design  freedom  and  have  been  analyzed  in  [5].  We  enploy  a  fiiU-wave  electromagnetic 
solver  based  on  this  formulation  to  predict  the  performance  of  each  shape.  In  general,  it 
is  more  difficult  to  design  an  effective  absorber  for  the  horizontal  polarization  under 
oblique  incidence  than  for  the  vertical  polarization.  Here  the  horizontal  polarization  is 
considered  using  the  standard  MAGRAM  coating  material  backed  by  a  conducting 
surfece.  We  compare  the  GA-optimized  shape  to  the  conventional  planar  and  triangular 
shaped  coatings.  Fig.  6  shows  the  corresponding  reflection  coefficients  over  a  frequency 
range  of  interest  for  the  three  different  shaped  coatings.  The  planar  shaped  coating  shows 


Fig.  5.  Geometry  of  the  corrugated  absorber  under  oblique  incidence. 
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Absorbing  Performance 


the  largest  reflection  (about  -5  dB)  within  the  frequency  band.  By  using  the  triangular 
shaped  profile  the  reflection  coeflScient  can  be  reduced  to  -lOdB.  The  GA-optrmized 
coating  shows  an  even  better  absorbing  performance  than  the  triangular  design.  The 
reflection  coeflScient  is  less  than  -15  dB  over  nearly  the  entire  frequency  range.  It  is  also 
noted  that  the  GA-optimized  shape  resembles  a  rectangular  profile.  One  possible 
interpretation  is  that  a  horizontally  polarized  wave  becomes  nearly  vertically  polarized  on 
the  vertical  sidewalls  of  the  rectangular  shaped  groove,  and  thus  is  easier  to  get  absorbed 
by  the  coating. 
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Fig.  7.  Pareto  front  showing  absorbing  performance  vs.  height  of  the 
coating  profile. 


Two  important  objectives  in  the  design  of  an  absorber  are  good  absorbing 
performance  and  low  coating  profile  (which  translates  into  light  weight).  Next,  we  utilize 
the  Pareto  GA  [6]  to  carry  out  this  multi-objective  optimization  eflBciently.  All  the 
samples  of  the  population  are  ranked  using  the  non-dominated  sorting  method  [7].  Based 
on  the  rank,  a  reproduction  process  is  performed  to  refine  the  population  into  the  next 
generation.  In  order  to  avoid  the  solutions  from  converging  to  a  single  point,  we  perform 
a  sharing  scheme  described  in  [8]  to  generate  a  well-dispersed  population. 

We  have  applied  the  Pareto  GA  to  corrugated  coating  design.  In  Fig.  7,  the  Pareto 
fi'ont  is  plotted  in  terms  of  the  absorbing  performance  versus  the  height  of  the  profile. 
Also  shown  in  the  insets  are  four  optimized  coating  shapes  that  have  different  heights. 
Inset  shape  (a)  shows  the  lowest  profile  among  the  four  samples,  but  has  the  highest 
reflection.  Inset  shape  (d)  has  the  highest  profile  and  the  lowest  reflection  among  the 
four  samples.  This  result  forms  a  useful  design  chart  for  signature  control  engineers  in 
trading  off  absorber  performance  against  coating  height. 

Electrically  Small  Wire  Antenna  Design.  The  design  of  electrically  small  antennas 
is  currently  of  great  interest  for  both  HF  communications  (where  antennas  are  on  the 
order  of  tens  of  meters)  and  wireless  handheld  devices.  However,  miniaturization 
impacts  both  antenna  efficiency  and  bandwidth,  two  critical  parameters  in  high  data  rate, 
low  power  consumption  systems.  Fig.  8  plots  the  well-known  theoretical  bandwidth  limit 
for  miniaturized  antennas  derived  by  L.  J.  Chu  in  1948  [9,10].  It  shows  that  as  the 
electrical  size  of  the  antenna  (measured  in  terms  of  hr  where  r  is  the  radius  of  the  smallest 
sphere  enclosing  the  antenna  and  k= 2 ^(wavelength))  is  decreased,  the  bandwidth  of  the 
antenna  is  also  reduced.  A  similar  trend  is  also  expected  for  the  efficiency  of  the  antenna, 
as  the  volume  over  which  the  current  on  the  antenna  can  flow  is  decreased.  Altshuler 
first  reported  on  the  use  of  GA  for  designing  electrically  small  wire  antennas  [11].  We 
have  applied  GA  in  conjunction  with  the  Numerical  Electromagnetics  Code  to  design 
electrically  small,  shaped  vdre  antennas  having  maximum  bandwidth  and  efficiency.  Our 
objective  is  to  devise  design  methodologies  that  can  approach  the  fundamental  limit  in 
miniaturization.  In  our  initial  design,  we  used  a  single-arm  wire  with  a  number  of  bent 
segments.  Fig.  8  shows  the  resulting  bandwidth  performance  fi-om  our  design  method 


Fig.  9.  Bandwidth  and  efficiency  of  the  kr=0.5  GA  antenna.  The  measured  half-power 
bandwidth  is  8.8%  and  the  measured  efficiency  (using  the  Wheeler  cap  method)  is  94%  at 
the  operating  frequency  of 400  MHz. 


versus  the  Chu  limit.  Fig.  9  shows  the  bandwidth  and  efiBciency  comparison  between  the 
designed  and  measured  results  of  a  particular  single-arm  GA  design  (^r=0.5).  While 
these  initial  design  results  are  encouraging,  they  are  still  significantly  below  the  Chu 
limit.  We  believe  additional  innovations  can  be  devised  to  better  approach  the  Chu  limit. 
Furthermore,  the  efiBciency  of  the  antenna  must  also  be  considered  in  the  miniaturization, 
and  the  Pareto  GA  scheme,  which  we  have  utilized  effectively  in  the  corrugated  absorber 
problem,  should  be  well  suited  to  explore  this  multi-objective  problem. 

C.  FOLLOW-UP  STATEMENT: 

We  have  obtained  very  encouraging  results  in  a  number  of  research  topics.  In  the 
coming  year,  we  will  continue  our  research  along  the  topics  outlined  above.  In  particular, 
we  will:  (i)  extend  our  work  in  microstrip  antennas  to  investigate  antenna  miniaturization 
using  patch  shapes  and  shorting  pins,  (ii)  investigate  the  effect  of  mutual  coupling  in 
microstrip  arrays  and  find  ways  to  mitigate  coupling  via  patch  shapes  to  improve  scan 
performance,  (iii)  fiirther  improve  the  design  methodology  of  electrically  small  wire 
antennas  to  approach  the  Chu  limit  and  map  out  performance  bounds,  (iv)  investigate 
platform  effects  in  smaU  antenna  design,  (v)  apply  Pareto  GA  to  more  fully  explore 
multi-objective  antenna  optimization  problems,  (vi)  e5q)lore  the  design  of  planar,  ultra- 
wideband  antennas. 

While  GA  is  an  efiBcient  global  optimizer  ideally  suited  for  exploring  very  complex 
design  problems  such  as  those  outlined  above,  it  typically  requires  a  very  large  number  of 
cost  function  evaluations.  We  will  leverage  against  our  previous  research  to  minimize 
the  number  of  CEM  calculations  needed  for  cost  function  computation.  Particular 
attention  will  be  paid  to  the  close  coupling  between  GA  and  the  CEM  simulator  to 
establish  a  framework  in  which  the  design  process  can  be  streamlined.  The  potential 
impact  of  the  proposed  research  is  that  novel  design  concepts  will  be  uncovered  that  can 
significantly  outperform  conventional  designs  for  sh^board  antenna  problems. 
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nisms  arising  from  multiple  scattering  (see  Fig.  1).  The  time 
of  arrival  and  the  strength  of  the  different  mechanisms  are 
extracted  from  the  available  frequency  data  using  the  adap¬ 
tive  feature  extraction  (AFE)  algorithm  [2].  It  was  shown  that 
good  interpolation  results  can  be  obtained  even  for  very 
sparsely  sampled  data  in  frequency. 

In  this  paper,  we  set  out  to  improve  the  accuracy  of  the 
above  interpolation  algorithm  for  complex  targets  by  incorpo¬ 
rating  additional  scattering  physics  in  the  model.  The  multi¬ 
ple-arrival  model  is  well  matched  to  the  ray-optical  descrip¬ 
tion  of  fields  and  currents.  However,  more  general  targets 
contain  not  only  large-scale  features,  but  also  small  resonant 
features.  Although  these  resonant  contributions  are,  in  gen¬ 
eral,  weaker  than  the  returns  from  large-scale  features,  they 
do  impact  the  accuracy  of  the  interpolation  if  not  properly 
taken  into  account.  The  multiple-arrival  model,  however,  is 
not  a  good  model  for  resonance  phenomena  as  many  time- 
of-arrival  terms  are  required  to  adequately  describe  reso¬ 
nance.  A  much  better  model  is  the  rational-function  model, 
which  has  been  widely  used  in  frequency-interpolation  prob¬ 
lems  in  the  resonant  region  [3-8],  It  can  be  shown  that  the 
multiple-arrival  model  and  the  rational-function  model  are, 
in  fact,  complementary  models.  The  former  is  a  sum  of 
exponentials  in  frequency,  while  the  latter  is  a  sum  of  expo¬ 
nentials  in  time.  To  take  advantage  of  the  strength  of  each 
model,  we  propose  here  a  hybrid  interpolation  scheme  that 
uses  the  most  appropriate  model  over  different  portions  of 
the  target.  At  each  point  on  the  target,  we  parameterize  the 
induced  current  using  both  models  separately.  Next,  by  using 
the  resulting  parameterization  error  as  the  selection  crite¬ 
rion,  we  choose  the  model  that  best  matches  the  scattering 
physics  at  that  location.  The  total  scattered  field  is  then 
constructed  from  the  interpolated  currents.  Section  2  de¬ 
scribes  our  formulation.  The  two  models  and  the  methods  for 
parameterization  are  first  summarized.  They  are  followed  by 
the  hybrid  procedure.  Our  test  results  in  Section  3  show  that 
the  hybrid  scheme  is  superior  to  either  the  multiple-arrival  or 
the  rational-function  model  alone  for  a  target  containing 
complex  features. 


2.  FORMULATION 

2.1,  Multiple-Arrival  Model  and  AFE  Interpolation  Algorithm. 
In  the  multiple-arrival  model,  we  assume  that  the  induced 
current  at  a  point  on  the  target  can  be  written  as  a  summa¬ 
tion  of  multiple-scattering  mechanisms  from  M  different 
incident  paths  (see  Fig.  1): 

/(/)=  EAexp(-;27r/t,)  (1) 

where  /  is  the  frequency,  t^  is  the  time  of  arrival  for  mecha¬ 
nism  /,  and  Aj  is  the  corresponding  excitation  amplitude. 
This  time-of-arrival  model  is  based  on  high-frequency  ray- 
optical  phenomena,  and  has  been  well  utilized  in  the  electro¬ 
magnetics  community  [9~11].  The  last  term  in  this  model 
incorporates  an  additional  frequency-dependent  factor  a,, 
which  was  not  used  in  [1].  It  is  consistent  with  high-frequency 
diffiraction  theory  [10,  12,  13],  and  improves  the  accuracy  of 
the  model.  Lastly,  the  extra  factor  in  the  denominator  is 
included  for  normalization. 

To  determine  the  unknowns  in  the  model  based  on  the 
available  frequency  samples  of  the  current,  we  use  the  adap¬ 
tive  feature  extraction  algorithm  [2].  The  parameterization  is 
carried  out  in  an  iterative  manner  starting  from  the  strongest 
mechanism.  The  sampled  frequency  response  is  first  pro¬ 
jected  onto  the  complex  conjugate  of  the  model  bases  for  all 
possible  values  of  tj  and  aj.  We  then  select  the  basis  that 
gives  the  maximum  projection  value.  This  is  described  as 
follows: 


Ai 


=  max( 


//(/),  exp(y27r^p 


(2) 


where  the  inner  product  in  the  above  formula  is  defined  as 

<^(/),M/))  =  -^Efl(/pM/p.  (3) 

After  the  strongest  feature  with  (A^,  a,)  is  captured,  that 

feature  is  subtracted  from  the  signal  to  generate  a  remainder 
signal: 


i(/)  =  MP  -  exp( (4) 

The  above  procedure  is  iterated  to  extract  the  parameters  for 
each  scattering  mechanism  until  the  remaining  signal  reaches 
a  sufficiendy  small  level. 

A  key  concept  in  AFE  is  to  use  random  frequency  sam¬ 
pling  during  the  original  data  generation  to  avoid  the  ambigu¬ 
ity  in  selecting  the  strongest  feature  [1].  Therefore,  the  fre¬ 
quencies  at  which  the  electromagnetic  computations  are 
carried  out  should  not  be  evenly  spaced.  In  addition,  since 
the  strongest  feature  contains  interference  from  the  weaker 
features,  the  amplitude  of  each  feature  is  not  extracted 
perfectly  using  ATO.  We  have  found  that,  by  performing  a 
linear  least  square  optimization  for  the  amplitude  parameter 
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after  the  whole  iteration  procedure,  the  AFE  performance  is 
significantly  improved. 


2.2.  Rational  Function  Model  and  Cauchy  Method  The  multi¬ 
ple-arrival  model  and  AFE  can  be  used  to  effectively  carry 
out  frequency  interpolation  for  targets  containing  ray-optical 
phenomena.  For  target  resonance,  however,  it  is  more  appro¬ 
priate  to  use  a  rational-function  model  to  describe  the  fre¬ 
quency  resonance  mechanism.  The  rational-function  model  is 
introduced  below,  and  the  Cauchy  method  for  determining 
the  model  parameters  is  summarized.  The  reader  is  referred 
to  [5,  6]  for  a  detailed  description  of  the  method. 

Consider  the  frequency  response  of  a  system  His).  From 
circuit  theory,  it  is  customary  to  describe  His)  by  the  ratio  of 
two  poljmomials  Ais)  and  Bis)  as  follows: 


His)^ 


k~0 

Q 


k~0 


(5) 


For  Q>  P,  this  is  a  pole  model,  and  thus  can  be  used  to 
describe  the  strong  resonance  behavior  of  the  current  in  the 
frequency  domain  by  the  proper  choice  of  the  pole  locations. 

Given  the  value  of  His)  and  its  frequency  derivatives  at 
some  frequency  points  s„,  the  Cauchy  problem  is  stated  as 

given  HKs„),  for  ;  =  0, . . . ,  n  =  1, . . . , 
find  P  and  Q,  =  0,..., P)  and  k  =  0, 

where  HKs„)  represents  the  yth  firequency  derivative  of  H  at 
Equation  (5)  can  be  rearranged  into  a  matrix  form  as 
follows: 


[A 


(6) 


where 


Aiuyk) 


k\ 

ik  ~i)! 


sll'^^^uik  -j) 


(7) 


and 


Here,  uik)  is  0  for  A:  <  0  and  1  otherwise,  Cjj  =  (yV(^-0  ““ 
0!)),  a  =  «/>F,and  b  -  ^of.Tosolve 

for  the  polynomial  coefficients  and  bj^  from  the  matrix 
equation  (5),  singular  value  decomposition  is  utilized.  The 
model  orders  P  and  Q  can  be  estimated  from  the  number  of 
significant  singular  values.  The  singular  vector  corresponding 
to  the  smallest  singular  value  is  chosen  to  be  the  solution  of 
(6). 

Several  comments  are  in  order.  First,  the  quantity  to  be 
interpolated  in  our  problem  is  the  frequency  response  of  the 
current.  Thus,  His)  ==  /(/)  where  s  =  y27r/.  Second,  we  only 
use  frequency  points,  and  not  frequency  derivatives  as  the 
input  to  the  Cauchy  algorithm.  This  is  for  ease  of  hybridiza¬ 
tion  with  the  AFE  algorithm.  Furthermore,  for  an  iterative 
solution  to  electromagnetic  integral  equations,  the  computa¬ 
tional  complexity  to  generate  the  frequency  derivative  of  the 


current  is  as  high  as  that  for  a  new  frequency  point.  Third, 
instead  of  estimating  the  orders  of  the  polynomial  functions, 
we  choose  g  =  P  +  1  and  P-l-Q-l-l=iVin  our  implemen¬ 
tation.  We  find  that  our  results  are  not  very  sensitive  to  the 
order  estimation.  Lastly,  we  note  that  the  frequency  points 
need  not  be  equally  spaced  in  carrying  out  the  Cauchy 
method.  This  gives  us  the  possibility  to  use  the  same  set  of 
sampled  data  for  both  the  Cauchy  and  AFE  interpolation 
schemes. 

2.1  Hybrid  AFE-Cauchy  Algorithm.  We  are  now  armed  with 
two  interpolation  algorithms.  The  AFE  algorithm  works  well 
for  low-order  multiple-scattering  mechanisms  that  resemble 
high-frequency  ray  optics,  while  the  Cauchy  algorithm  is 
better  suited  for  describing  resonance  phenomena.  The  es¬ 
sential  idea  of  the  hybrid  algorithm  is  to  choose  the  most 
appropriate  model  for  the  current  at  each  point  on  the  target. 
Our  approach  is  as  foUows.  At  each  point  on  the  target,  we 
first  parameterize  the  induced  current  using  both  the  AFE 
and  Cauchy  algorithms  separately,  with  the  same  set  of 
sampled  frequency  responses.  The  model  that  best  matches 
the  scattering  physics  at  that  location  is  then  chosen.  In  the 
actual  implementation,  the  selection  is  made  by  comparing 
the  interpolated  data  from  the  two  models  against  the  brute- 
force  computation  at  several  new  frequency  points.  The  model 
with  the  smaller  interpolation  error  is  chosen  for  that  loca¬ 
tion.  As  a  result,  the  surface  of  the  target  can  be  divided  into 
two  regions — the  AFE  region  and  the  Cauchy  region.  Over 
the  AFE  region,  the  currents  exhibit  high-frequency  scatter¬ 
ing  phenomenology.  Over  the  Cauchy  region,  the  currents 
exhibit  strong  resonance.  Once  the  parameterization  and 
model  determination  have  been  made,  the  interpolated  cur¬ 
rents  on  the  target  over  the  frequency  band  of  interest  can  be 
generated.  Therefore,  the  scattered  far  field  can  be  obtained 
by  integrating  the  induced  current  over  the  target  at  any 
frequency  of  interest. 

It  has  been  pointed  out  that  the  AFE  algorithm  requires 
as  its  input  the  currents  computed  at  a  set  of  nonuniformly 
sampled  frequencies.  In  fact,  the  more  random  the  sampling, 
the  better  the  performance  of  the  AFE.  It  has  also  been 
noted  that,  in  the  Cauchy  algorithm,  the  frequency  points 
need  not  be  uniformly  sampled.  However,  equally  spaced 
frequency  points  increase  the  chance  for  the  Cauchy  method 
to  accurately  capture  the  poles  in  the  frequency  domain. 
Consequently,  we  find  contradictory  requirements  for  the 
optimal  sampling  in  frequency  between  the  AFE  and  Cauchy 
algorithms.  AFE  requires  random  sampling,  while  Cauchy 
prefers  uniform  sampling.  Here,  we  choose  a  compromise 
solution,  i.e.,  semirandom  sampling.  To  choose  the  N  fre¬ 
quency  points  needed  for  interpolation,  we  first  equally  divide 
the  whole  frequency  band  into  N  uniformly  spaced  subbands. 
The  sampling  point  within  each  subband  is  chosen  randomly 
based  on  a  uniform  probability  distribution.  The  sampled 
data  can  then  be  used  in  both  algorithms. 

3.  NUMERICAL  EXAMPLE 

We  consider  a  2-D  conducting  target  as  shown  in  Figure  2. 
Although  this  plate-like  target  looks  simple,  its  scattering 
characteristics  involve  a  variety  of  scattering  mechanisms  [14]. 
An  £-polarized  plane  wave  is  incident  at  an  angle  of  60®  from 
the  right.  The  reference  backscattering  data  versus  frequency 
are  first  computed  using  the  method  of  moments  (MoM)  at 
100  frequency  points  from  0.1  to  10  GHz  in  0.1  GHz  steps. 
The  frequency  response  of  the  backscattered  field  is  plotted 
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Rgure  2  Geometry  of  the  conducting  plate  with  cavity  and  fin 


as  the  solid  line  in  Figure  3.  By  inverse  Fourier  transforming 
the  frequency  data,  the  range  profile  of  the  target  can  be 
generated.  It  is  plotted  in  Figure  4  as  the  solid  curve  (a 
Hamming  window  has  been  applied  to  reduce  range  side- 
lobes).  We  recognize  four  large  peaks  in  the  range  profile, 
which  correspond,  respectively,  to  the  scattering  from  the 
right  edge  of  the  plate,  the  aperture  of  the  partially  open 
cavity,  the  interior  of  the  cavity,  and  the  small  fin  at  the  left 
edge.  The  resonant  behavior  of  the  cavity  interior  can  clearly 


ie"i 
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Frequency  (GHz) 
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Figure  3  (a)  Comparison  of  the  backscattered  field  versus  fre¬ 
quency  between  the  brute-force  MoM  and  the  AFE  interpolation 
results  based  on  20  points,  (b)  d^mparison  of  the  backscattered  field 
versus  frequency  between  the  brute-force  MoM  and  the  Cauchy 
interpolation  results  based  on  20  points,  (c)  Comparison  of  the 
backscattered  field  versus  frequency  between  the  brute-force  MoM 
and  the  hybrid  interpolation  results  based  on  20  points 


Range 

(c) 


Rgure  4  (a)  Comparison  of  the  range  profiles  generated  from  the 
brute-force  MoM  and  the  AFE-interpolated  data,  (b)  (Ilomparison  of 
the  range  profiles  generated  from  the  brute-force  MoM  and  the 
Cauchy-interpolated  data,  (c)  Comparison  of  the  range  profiles  gen¬ 
erated  from  the  brute-force  MoM  and  the  hybrid-interpolated  data 


be  seen  as  dispersive  ringing  occurs  after  the  third  peak  and 
overlaps  with  the  fourth  peak.  We  should  note  that  the 
frequency  response  is  marginally  sampled  at  a  step  size  of  0.1 
GHz.  Further  undersampling  will  result  in  a  loss  of  features 
in  the  frequency  domain  or  aliasing  in  the  time/range 
domain. 

We  next  cany  out  the  frequency  interpolation  using  only 
20  frequency  points.  The  locations  of  the  20  frequency  points 
are  chosen  semirandomly,  as  discussed  in  the  last  section,  and 
are  indicated  as  small  arrows  along  the  frequency  axis  in 
Figure  3.  First,  the  target  current  is  interpolated  using  the 
AFE  algorithm  based  on  the  multiple-arrival  model  over  the 
whole  target  surface.  The  interpolated  result  is  plotted  as 
the  dashed  line  in  Figure  3(a).  Compared  to  the  reference 
result  from  the  brute-force  MoM,  AFE  gives  rather  good 
performance,  and  can  capture  most  of  the  features  in  the 
scattered  field.  The  largest  error  in  the  AFE  result  occurs  at 
frequencies  around  4.5  and  9  GHz,  which  are  close  to  the 
first  two  resonant  frequencies  of  the  small  cavity.  This  obser¬ 
vation  is  further  corroborated  in  the  corresponding  range 
profile  from  the  AFE-interpolated  data,  shown  as  the  dashed 
curve  in  Figure  4(a).  The  AFE  algorithm  performs  well  in 
predicting  the  three  strong  peaks  due  to  the  exterior  features 
on  the  target.  However,  it  has  trouble  predicting  the  resonant 
cavity  return,  which  is  rather  dispersive  in  range. 
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Next,  the  Cauchy  method  is  used  to  interpolate  the  cur¬ 
rent  over  the  entire  target  surface.  The  frequency  response 
and  the  range  profile  results  are  shown,  respectively,  in 
Figures  3(b)  and  4(b)  as  dashed  lines.  Since  the  rational 
function  is  a  pole  model,  we  find  that  the  Cauchy  algorithm 
tries  hard  to  capture  the  poles  of  the  system.  However,  since 
the  current  is  sparsely  sampled  in  frequency,  an  incorrect 
pole  location  and  pole  strength  are  generated.  This  is  particu¬ 
larly  true  for  low-order  multiple-scattering  events  that  have  a 
large  delay  spread  in  the  times  of  arrival.  The  long  delay 
spread  in  time  translates  into  rapid  oscillations  in  the  fre¬ 
quency  response,  and  a  very  high  model  order  is  needed  to 
adequately  model  the  response  using  poles.  From  the  range 
profile  in  Figure  4(b),  we  see  that  the  Cauchy  algorithm  gives 
extremely  poor  results.  Thus,  the  Cauchy  algorithm  cannot  be 
used  by  itself  to  carry  out  the  frequency  interpolation  of  large 
targets  that  are  dominated  by  ray-optical  phenomena. 

Finally,  we  apply  the  hybrid  interpolation  algorithm.  The 
selection  between  the  Cauchy-interpolated  and  the  AFE-in- 
terpolated  currents  at  each  point  on  the  target  surface  is 
made  by  comparing  them  against  the  reference  current  at 
three  extra  test  frequencies.  The  chosen  interpolated  cur¬ 
rents  are  then  used  to  generate  the  backscattered  field,  as 
shown  in  Figure  3(c).  When  compared  to  Figure  3(a)  and  (b), 
the  hybrid  result  shows  an  obvious  improvement  in  its  agree¬ 
ment  with  the  reference  result.  Note  that  the  error  near  the 
cavity  resonant  frequencies  in  Figure  3(a)  is  greatly  reduced 
in  Figure  3(c).  The  range  profile  in  Figure  4(c)  further  sub¬ 
stantiates  this  claim  as  both  the  scattering  centers  and  the 
resonance  of  the  cavity  are  well  predicted. 

Figure  5  shows  how  the  whole  target  is  automatically 
divided  into  two  regions  during  the  hybrid  procedure,  de¬ 
pending  on  the  local  scattering  physics.  The  Cauchy  region  is 
shown  within  the  open  lines,  while  the  AFE  region  is  shown 
within  the  solid  lines.  The  right  half  of  the  plate  fits  the 
multiple-arrival  model  very  well,  and  thus  is  chosen  as  the 
AFE  region.  The  cavity  in  the  middle  exhibits  strong  reso¬ 
nance,  and  can  be  well  parameterized  by  the  Cauchy  method. 
'^As  for  the  left  half  of  the  plate,  we  expect  it  to  be  an  AFE 
region  with  the  same  scattering  mechanisms  as  the  right-half 
part.  However,  since  the  delay  spread  in  the  times  of  arrival 
between  the  direct  incident  wave  and  the  scattered  wave  from 
the  fin  at  the  left  edge  is  small,  the  Cauchy  method  has  a 


Cauchy  region 
AFE  region 


Figure  5  Illustration  of  how  different  portions  of  the  target  are 
parameterized  in  the  hybrid  approach.  The  Cauchy  region  is  shown 
within  the  open  lines,  and  the  AFE  region  is  shown  within  the  solid 
lines 


sufficient  model  order  to  parameterize  the  scattering  behav¬ 
ior.  In  fact,  we  find  that  either  the  Cauchy  or  the  AFE  model 
may  be  chosen  in  this  region  with  equally  good  results. 

4.  SUMMARY 

In  this  paper,  we  have  proposed  a  frequency-interpolation 
algorithm  for  electromagnetic  scattering  data  from  large, 
complex  targets.  It  is  based  on  the  hybridization  of  the 
existing  multiple-arrival  model  and  the  rational-function 
model.  The  adaptive  feature  extraction  algorithm  is  used  to 
extract  the  parameters  in  the  multiple-arrival  model,  while 
the  parameters  in  the  rational-function  model  are  obtained 
using  the  Cauchy  method.  At  each  point  on  the  target,  we 
first  parameterize  the  induced  current  using  both  the  AFE 
and  Cauchy  algorithms  separately.  The  model  that  best 
matches  the  scattering  physics  at  that  location  is  automati¬ 
cally  chosen  based  on  the  resulting  interpolation  error.  Con¬ 
sequently,  those  regions  on  the  target  that  are  best  described 
by  ray-optical  phenomena  are  interpolated  by  the  multiple- 
arrival  model,  while  those  regions  that  exhibit  resonance 
phenomena  are  interpolated  by  the  rational-function  model. 
Numerical  results  for  a  target  containing  complex  features 
have  demonstrated  that  the  hybrid  model  results  in  more 
accurate  interpolation  than  either  of  the  two  models  alone. 
The  hybrid  interpolation  scheme  is  quite  robust,  and  can  lead 
to  significant  computational  savings  since  broadband  signa¬ 
ture  data  can  be  generated  from  a  very  sparse  set  of  com¬ 
puted  frequency  points.  Extension  of  this  algorithm  to  3-D, 
nonconducting  targets  is  currently  being  investigated. 
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Paper  [13] 


Wavelet-lMised  preconditioner  for  threo- 
dimensional  electromagnetic  integral 
equations 

R  Deng  and  H.  Ling 

A  wavcict*based  nKthod  is  piopo^  to  cflbctivc^  pfcooncfitton 
3D  doctromagaedc  integipBl  eq[iialioii$.  The  appmximaicHhnmase 
precopdMoner  a  oonstracted  in  the  wavdet  dotnoin  where  both 
the  tmuienl  nsOm  and  its  mveoD  cxli%a  ipoise; 
sltiiciuie&  The  invcnioii  is  caniod  out  ai  a  Frobenhis^ioim 
mioaniuitifio  {woblem.  Numcricai  icsulU  on  a  3D  cavity  show 
that  the  itcnlkn  numbos  aie  i^gnificanify  letfiiODd  with  the 
preoonefitioaed  system.  The  oompmalional  cost  of  the 
preconditifMicrg  hq)t  under  OCMogW). 

bttrodmtUtm:  Tbeie  ts  grawmg  inlctcsl  in  Ibe  con^xitalional  doo- 
troraeoBtks  oniHioiiily  im  the  iw  of  the  fiis^ 

OFMM)  [1]  for  Ecltmg  bfse-scab  decUnmagticlic  int^mol  eqiMr 
tions;  With  the  nmldiBvd  wiplementatinn  of  the  FMM,  the  oom- 
putotinnrf  comptouly  and  atorage  leqnhenicut  can  be  icAiced  to 
C<iVlog?Q  fix  jte  (moment  imlinHwtoid  product  To  take 
advantage  of  this  mfaned  complraaty,  iterative  sotvera  miKt  be 
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used.  However,  the  convergence  rate  of  the  iterative  solution  proc¬ 
ess  is  strongly  problera-depwident.  For  scattering  geometries 
involving  multiple  interactions  such  as  partially  open  cavities,  the 
moment  matrices  are  severely  ill-conditioned  and  preconditioning 
is  need  to  accelerate  the  convergence  rate  [2].  Some  well-docu¬ 
mented  preconditioning  methods  such  as  incomplete  LU-factorisa- 
tion  (ILU)  [3]  may  be  cfTectivc,  but  usually  require  wdl-above 
0(MogAf)  complexity  to  implement  In  this  Letter,  we  propose  a 
new  wavelet-based  method  to  construct  an  effective  preconditioner 
for  3D  moment  equations.  The  algorithm  is  based  on  determining 
an  approximate  inverse  for  the  moment  matrix  in  the  wavelet 
transform  domain.  The  computational  cost  of  the  preconditioncr 
is  kept  under  0(Mog^,  making  it  compatible  with  the  FMM. 

Problem  formulation:  Consider  a  3D  moment  equation  in  the 
space  domain: 

[Z]J  =  E  (1) 

where  [Z],  J  and  E  are,  respectively,  the  momait  matrix,  the 
induced  current  vector  and  the  excitation  vector.  If  eqn.  I  is  left- 
preconditioned  by  a  preconditioncr  [P],  we  obtain 

[P]IZ]J  =  [P]B  (2) 

Our  approadi  is  to  try  to  construct  a  [P]  that  is  the  approximate 
inverse  of  [2^  uang  the  wavelet  Iransfonn.  The  moment  equation 
in  eqn.  1  can  be  represented  using  wavelet  ba^  as 

[Z]J  =  E  (3) 

where 

[Z]  =  [M]''’12][M1  [Jl  =  [Ml’-J  m  = 

(4) 

and  [IVf]  is  the  unitary  wavelet  transform  matrix.  Note  that  in  30 
problems  [Z)  is  strongly  diagonai-dominant  Cons^uentiy, 
can  be  effectively  approjinated  by  a  sparse  matrix  with  the  multi¬ 
level  Tingeri  pattern  discussed  in  [4].  Furthermore,  we  observe 
from  eqn.  4  that 

[Z]”‘  =  [M]"’[ZrMMl  =  lZ]"'  (5) 

Since  the  smooth  parts  in  [Zf*  arc  converted  into  small  wavelet 
coefficients  through  the  transform,  [Z]-‘  can  again  be  approx- 
mated  by  a  sparse  matrix  with  the  same  ffn^r  pattern  as  {Z]. 
With  the  chosen  sparsity  pattern  for  teth  [Z]  and  [Z] the 
approximat&'inverse  preconditioner  for  [ZJ  can  now  be  solved  effi¬ 
ciently  by  castmg  it  as  a  Frobenius  norm  minimisation  problem 

[5]: 

<«) 

where  1  is  the  identity  matrix.  Since 

|lPllZ]-lL  =  E||Pi[Z)-eif  (7) 

i=i 

where  and  ^  are,  respectively,  the  fih  rows  of  [P]  and  I,  the 
solution  of  eqn.  6  beoennes  the  following  N  independent  least 
squares  problems 

min||Pj[Z]-e,|p  j  =  l,2,..„JV  (8) 

Once  [P]  is  solved  for  in  the  wavelet  domain,  the  approximate- 
inverse  preconditioner  [P]  can  be  obtained  by  the  inverse  wavelet 
transform: 

lP]  =  [M][Pl[M]^=.[Zr^  (9) 

Combining  eqns.  2  and  9  vre  have  the  following  preconditioned 
moment  equation: 

[M][P][M]^[Z1J  =  [M][P][M]^E  (10) 

To  summarise,  the  wavelet  preconditioncr  [P]  is  generated  by  first 
finding  the  approximate  wavelet  transform  of  the  moment  matrix 
IZ]  using  eqn.  4,  and  then  soK^  the  Frobenius-norm  minimisa¬ 
tion  problem  in  eqn.  8.  Once  [P]  is  obtained,  the  preconditioning 
operation  is  carried  out  by  a  scries  of  matrix-vector  multiplication 


operations  in  eqn.  10.  We  now  consider  the  computational  cost  of 
constructing  [P]  and  that  for  the  preconditioning  operation.  To 
construct  [P],  we  must  obtain  an  approximate  [Z],  Owing  to  the 
strong  source  sii^ularity  in  3D  problems,  [Z|  can  be  made  approx¬ 
imately  sparse  with  0(N)  nmizero  elements  by  ignoring  the  far- 
field  interactions.  Therefore,  the  conqnitation  of  can  be  imple¬ 
mented  with  fast  t^/^haimcl  filterbank  filtering  with  about  0(N) 
operations.  Once  [Z]  is  obtained,  the  complexity  to  solve  eqn.  8 
can  be  made  prop^onai  to  the  |MT>blcm  aze  JV,  ^ovided  that  we 
carefully  control  the  sparsity  patterns  of  [Z]  and  its  inverse.  Next 
we  consider  the  complexity  to  carry  out  the  preconditioning  oper¬ 
ation  on  the  lell-hand  side  of  eqn.  10.  Since  there  are  only 
0(iVlog/y)  nonzero  elements  in  the  wavelet  transform  matrix  [Ml. 
the  product  of  [IVfl  with  any  vector  rcqtiires  0(N\ogN)  operations. 
Note  also  that  [P]  is  a  sparse  matrix  with  0{N)  nonzero  elements. 
Therefore,  if  the  preconditioning  operation  is  carried  out  as  a 
series  of  matrix-vector  products  on  [Z]J  or  E,  the  total  computa¬ 
tion  cost  is  limited  to  0(Mo$N)  per  iteration. 


Logarithmic  scale  with  L  =  4  (IV  =  1024) 
a  Wavelet  transformed  moment  matrix 
b  Inverse  of  transformed  matrix 


Fig,  2  Non-zero  pattern  used  to  solve  Frobenius-norm  minimisation 
problem 

Numerical  results:  The  algorithm  is  tested  using  a  3D  conducting 
rectangular  cavity  with  an  open  end.  The  (width)  x  (height)  x 
(depth)  dimensions  of  the  cavity  me  Rx  Rx  1.5i{  whm  Risa 
size  parameter.  The  incident  pl^  wave  is  horizxmtally  polarised 
with  a  frequent^  of  5.9GHz  and  makes  an  an^e  of  45^  fixMn  the 
cavity  opening  in  the  devation  plane.  The  problem  is  formulated 
in  terms  of  an  dectrio-field  int^ral  equation.  The  Rao-Wilton- 
Glisson  basis  functions  are  used  to  form  the  origmal  spao&do>i)^ 
moment  equation  with  a  discretisation  size  of  about  VlO  (or 
0.2").  The  greyscale  magnitude  plots  of  the  wavdet-based  moment 
matrix  [Z]  and  its  inverse  [Z]"*  are  shown,  respectively,  in 
RgS.  la  and  b  with  a  problem  size  of  =  1024  (/I  «  2.1").  The 
wavelet  filter  used  in  the  transform  is  the  Daubechies  filter  of 
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order  6,  and  the  maximum  wavelet  tranrform  level  is  L  =  4.  As 
expected,  both  matiioes  exhibit  the  multikvd  *fingcr’  structure. 
Accordingly  we  propose  to  use  a  fixed  ^jaise  matrix  ^ttem 
shown  in  Fig.  2  with  finger  width  D  for  both  [2^  and  in 
solving  the  Frobenius-nonn  minimisation  pioWenL  If  the  param^ 
ters  D  and  L  arc  kqpt  constant,  the  oomp^ty  to  solve  eqn.  10  is 
proportional  to  the  problem  size  N,  To  investigate  the  perfom- 
ance  of  the  preconditioner,  we  proportionally  increase  the  physical 
size  of  the  cavity  with  the  parameter  R  changing  from  1.2”  to  5.1" 
so  that  N  increases  from  256  to  8192,  while  the  discretisation 
interval  and  the  frequency  remain  unchanged.  The  resulting  itera¬ 
tion  numbers  in  solving  the  moment  equations  as  a  function  of  the 
problem  size  are  shown  in  Fig.  3.  The  iterative  solver  is  BICC- 
STAB  and  a  relative  residual  error  of  0.001  is  used  as  the  stopping 
criterion.  The  maximum  wavelet  transform  level  is  Z,  =  5  and  the 
finger  width  D  of  the  sparsity  pattern  is  16  for  iV  <  1024  and  32 
for  N  >  1024.  We  observe  that  with  preconditioning,  the  iteration 
numbers  are  very  small  and  grow  at  a  very  slow  rate  with  respect 
to  the  problem  size.  The  results  demonstrate  the  effectiveness  of 
the  new  wavelet-based  proconditioner  for  3D  moment  equations. 


Fig.  3  Iteration  number  against  problem  size  for  cavity  problem 


— without  preconditioning 
— + —  with  wavelet  preconditioner 

Conclusions:  We  have  proposed  a  wavelet-based  preconditioner 
for  3D  electromagnetic  integral  equations  in  this  work.  Numerical 
results  showed  the  precondilioner  is  very  effective  for  ill-condi¬ 
tioned  cavity  structures.  The  total  computational  cost  for  the  pre¬ 
conditioning  can  be  kept  under  O(MogA0-  The  new  algorithm  is 
compatible  with  fast  boundary-integral  algorithms  such  as  the 
multilevel  FMM. 
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Stiape<iplinffeMtion  of  broa<fliand  micros^ 
aitennas  using  genetic  atgorithm 


H.  Choo.  A.  Hutani,  L.C.  Trintinalia  and  H.  ling 

Ibe  algoiiiha  (GA)  b  used  to  dcagn  po^ 
iitirjRnitifip  ai^tennns  on  FR<4  suhstiaic  for  bfotufoond 


K  of  foe  GArOptiBBKd  dcatgg 


htavduafon:  it  is  weQ  known  that  standard  mionoslrip  pafdi 
«hgA  voy  nanov  baadwicfth.  Varfoiis  bioadbandiBg 
yiy^hndk  havo  proposed  to  date.  For  mstance,  adifing  pa^ 
fHtff-  |Brteheg,  iiM»g  thipir  air  mhsiiatE;  Slnckil^  patrfjCS  and  IIStDg 
slMK&ie  pod  for  reacth«  loaflfotK  aie  wdi  faram  tediini^ 
^ffWnuCiie  micioatrip  baodiridtii  [ll*  Reoen^,  JofaDson  and 
Raianft-Shnii  icported  on  the  i»b  of  the  genetic  algqiithm  (GA) 
to  sesBch  for  novel  patdi  dnpes  for  beoadband  opentdras  PI 
The  attiacthwpeas  of  GA  shape  t^kniuitinn  is  diat  impniived 

pqr6wTftnn<y».g!nti  hoiMliiMwdwittMMttncreasipgo^ 
vidoM  or  numAchakig  oooC  Tbqr  md  duck  afr 
caqtoed  niela&  pntdi  SOBS  vp  to  half  a  wavelci^ 
fo  tUs  Letfoi;  w  caanm  the  loe  oTGA  for 
appGcatioaa  la  ocwtiasr  to  dm  wwk  of  Jtiaaoo  Rdnto- 
SShob*  we  cn^lofy  FR4  as  foe  sdxsbatB  material,  since  il  is  foe 
inntf  oomnotebr  used  material  m  wnriess  derices,  to  addidon, 
-  fewer  gpcgnetricai  comttaato  aic  used  m  foe  GA  m  hope  of 

ofeteinnis  better  ^tfoaloptinnBi^WgiqMrl  a  fagfeddbaadwidfo 

iimwowHnwiit  fiam  cor  GA-opfonoed  nicroslrip  foape  oanqMied 
to  that  of  a  standard  square  mktostito 
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GA  optimisation:  GA  is  implemented  to  optimise  the  microstrip 
patch  shape  to  achieve  broad  bandwidth.  In  our  GA,  we  use  a 
twoKlimensional  (2D)  chromosome  to  encode  each  patch  shape 
into  a  binary  map  [3].  The  metallic  sub-patdies  are  represented  by 
ones  and  the  no-mctal  areas  arc  represented  by  zeros.  Since  it  is 
more  desirable  to  obtain  optimised  patch  shapes  that  arc  well-con¬ 
nected  from  the  manufacturing  point  of  view,  a  2D  median  filter  is 
applied  to  the  chromosomes  to  create  a  more  realisable  population 
at  each  generation  of  the  GA. 

To  evaluate  the  performance  of  eadi  patch  shape,  a  full-wave 
electromagnetic  patch  code  is  used  to  pi^ict  its  bandwidth  per¬ 
formance  [4].  The  formulation  of  the  code  is  similar  to  that 
described  in  [5]  and  is  based  on  the  solution  to  the  etectric  field 
integral  equation  with  the  periodic,  layered  medium  Green’s  func¬ 
tion  as  its  kernel.  Roof-top  basis  functions  are  used  to  expand  the 
unknown  current  on  the  metal  patch  and  fast  Fourier  transform  is 
used  to  accelerate  the  computation  of  the  matrix  elements. 
Because  of  the  assumed  periodkaty  of  the  patches  in  this  code,  we 
use  a  large  enough  period  to  simulate  a  single  patch. 

The  design  goal  is  to  broaden  the  bandwidth  of  a  microstrip 
antenna  with  a  centre  frequency  of  2GHz  by  changing  the  patch 
shape.  To  achieve  the  design  the  cost  function  is  defined  as 
the  average  of  those  5ii  values  that  exceed  -lOdB  (i.e.  VSWR  = 
2:1)  within  the  frequency  band  of  interest  The  target  frequency 
range  is  between  1,9  and  ZlGHz. 

Based  on  the  cost  function,  the  next  generation  is  created  by  a 
reproduction  process  that  involves  crossover,  mutation  and  2D 
median  filtering.  A  two-point  crossover  scheme  using  three  chro¬ 
mosomes  is  devised.  The  process  selects  three  diromosomes  as 
parents,  and  divides  each  chromosome  into  three  parts.  Intermin¬ 
gling  the  three  parent  dutunosomes  then  makes  three  child  diro¬ 
mosomes.  This  crossover  scheme  exhibits  a  mote  disruptive 
characteristic  for  regeneration  than  the  conventional  one-point  or 
two-point  crossover.  It  serves  to  counteract  the  median  filtering 
effect  and  is  found  to  result  in  better  convergence  rate.  The  repro¬ 
duction  process  is  iterated  until  the  cost  ftjnction  converges  to  a 
minimum  value. 


Fig.  1  Schematic  diagram  and  return  loss  of  GA-optimised  microstrip 
antenna 

a  GA-optimised  microstrip  antenna  using  16  x  16  resolution  within  72 
X  72mm  area 

Grey  pixels  are  metal  and  white  dot  shows  position  of  probe  feed 
b  Return  loss  of  antenna 
—  simulation 
□  meusurement 

Results:  Fig.  la  shows  the  shape  of  tlie  GA-optimised  microstrip, 
A  72  X  72mm  square  design  area  in  which  the  metallic  patch  can 
reside  is  discretised  into  a  16  x  16  grid  for  the  chromosome  defini¬ 
tion.  The  thickness  of  the  FR-4  substrate  (dielectric  constant  of 
4.3)  is  1.6mm.  In  the  GA-optimised  shape,  the  grey  pixels  are 
metal  and  the  white  pixels  have  no  metal.  The  white  dot  shows  the 
position  of  the  probe  feed.  To  experimentally  verify  the  GA 
design,  we  have  built  and  measured  such  a  microstrip  patch. 
Fig.  shows  the  return  loss  comparison  between  the  measure¬ 
ment  and  simulation  results.  The  solid  line  is  the  measuiement 
result  taken  from  an  HP8753C  network  analyser.  The  square  dots 
represent  the  simulation  result  Good  agrcemfait  can  be  observed 
b^ween  the  measurenient  and  simulation.  The  graph  shows  a 
baiKlwidth  of  6.16%  by  simulation  and  6.18%  by  measurement 
This  is  about  three  times  that  of  a  square  mterostrip  antenna  (36  x 
36mm),  which  has  a  bandwidth  of  1.98%.  Further  improvements 
in  the  bandwidth  can  be  obtained  from  the  GA  by  increasing  the 
grid  resolution  from  16  x  16  to  32  x  3Z  Figs.  2a  and  b  show, 
req)ectively,  the  GA-optimised  patdi  shape  and  the  bandwidth 
performanoe  in  the  higher  resolution  design.  The  bandwidth  is 


found  to  be  8.04%  by  simulation  and  8.10%  by  measurement  This 
is  about  four  times  that  of  a  square  microstrip  antenna. 

Finally,  the  operating  principle  of  the  GA-optimised  shape  is 
int«pret^.  It  is  clear  from  the  two  frequency  dips  in  Fig.  2b  that 
the  antenna  contains  two  operating  modes  that  are  very  closely 
spaced  in  frequency.  We  have  verified  the  two  modes  by  examin¬ 
ing  the  current  distributions  on  the  patch  at  1.99  and  Z07GHz.  In 
addition  to  the  duabnode  principle,  another  important  bandwidth 
enhancement  effect  is  achieved  through  the  rag^  edge  shape.  We 
have  found  that  when  the  patch  is  restricted  to  singlemode  opera¬ 
tion  (by  imposing  symmetry  constraints),  the  introduction  of  rag¬ 
ged  edges  in  the  GA-optimised  shape  can  enhance  the  bandwidth 
by  -'30%.  Therefore,  the  GA-optimised  design  combines  both  the 
duahnode  operation  and  ragged  edge  shape  to  adiieve  the  broad¬ 
est  bandwidth. 


Fig.  2  Schematic  diagram  and  return  loss  of  GA-optimixed  microstrip 
antenna 

a  GA-optimised  microstrip  antenna  using  32  x  32  resolution  within  72 
X  72mm  area 

Grey  pixels  are  metal  and  white  dot  shows  position  of  probe  feed 
b  Return  loss  of  antenna 
—  simulation 
□  measurement 

Conclusions:  Optimised  patch  shapes  for  microstrip  antennas  on 
thin  FR-4  substrate  have  been  investigated  using  the  genetic  algo¬ 
rithm.  The  optimised  shape  shows  a  fourfold  improvement  in 
bandwidth  when  contrasUxl  with  a  standard  square  microstrip 
antenna.  This  result  has  been  verified  by  laboratory  measuiement 
The  basic  operating  principle  of  the  optimised  shape  can  be 
explained  in  terms  of  a  combination  of  dualmode  operation  and 
ragged  edge  shape. 
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ABSTRACT^  77ir  txx^j&ig  matrix  is  a  samdard  way  to  vepieseau  the 
muiuai-couphng  ej^cf  betweat  ekutaus  m  an  anUaaa  array,  bt  this 
paper,  two  conmumfy  used  iigpproaches  to  ihsmmae  the  eot^Sng  matrix 
are  examined,  and  Aeir  Bmitatkms  ate  dacussed.  The  coapBng  matrix 
cleaned  pom  die  imnmmmmean-squateerwr  matching  of  Aeaeiice 
and  stand-okmeekoKUt  patterns  is  shown  to  be  mote  e^^^x  than  the 
munudringiedance  approach.  The  coiqJingiuaoix  modei  tr  tdso  ex- 
tended  for  mom  congdexaatemaa  arrays,  and  an  exampte  is  proaded  to 
dbtstnae  die  egeedveness  of  the  extended  nuML  O  2001  Join  WSey  & 
Sons,  loc.  Microwawe  Opt  Tedmol  Lett  28: 231-237, 2001. 

Kxywordsi  amty  muma!  cotqi&ig:  cotqdmg  matrix; 
actice  ekmaa  patterns 

L  INTRODUCTION 

It  IS  well  known  that  the  radiadon/recciving  patteni  of  an 
antenna  element  situated  in  an  array  environment  can  be 
quite  different  from  its  stand-alone  pattern  due  to  mutual¬ 
coupling  effects.  The  antenna  element  pattern  in  the  array 
environment  is  called  the  active  efement  pattern  in  the  an¬ 
tenna  communiqr  [IX  In  qiplicatkHis  such  as  directkMt  find¬ 
ing,  the  active  element  patterns  of  die  array  elements  are 
needed  to  cany  out  dhection-of-anivai  esthnattoii.  If  the 
stand-alone  element  patterns  are  osed  instead  of  the  active 
element  patterns,  significant  degradation  can  result  in  array 
performance  [2]. 

To  properly  account  for  the  mntnaTaHifding  ^fect  and 
for  simplicity,  it  is  customary  to  express  the  active  demmit 
patterns  as  the  product  of  a  couidhig  matrix  and  the  stand¬ 
alone  dement  patterns.  Altfaoi^  the  ooiqiling^iatrix  con¬ 
cept  is  widely  utilized  fiar  array  s^nal  processto^  two  very 
different  approaches  have  been  adc^ited  in  the  literature  to 
mterpret  and  determine  the  ooiqplii^  matrix.  The  first  ap¬ 
proach  conies  out  of  the  antenna  communi^.  Gnpta  and 
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Ksienski  [3]  first  derived  the  coupling  matrix  from  a  mi¬ 
crowave  network  point  of  view.  The  array  is  viewed  as  a 
multiport  network,  and  it  is  shown  that  the  coupling  matrix 
can  be  simply  related  to  the  generalized  impedance  matrix  of 
die  muhqimt  network.  Since  the  generalized  impedance  ma¬ 
trix  is  an  intrinsic  property  of  the  antenna  structure,  sudi  a 
derivation  of  the  conpliiig  matrix  is  quhe  attractive  ffom  the 
fundamental  electromagnetics  point  of  view.  This  approadi 
has  been  utilized  and  extended  in  a  number  of  antenna 
mutual-coupling  studies  [4-6]. 

The  second  approach  or^inated  from  the  array  rignal- 
procesring  community.  Friecflander  and  Wmss  [7]  first  pro¬ 
posed  a  method  for  direction  finding  based  on  the  coupling- 
matrix  modcL  In  thdr  wew,  fbc  oouf^g  matrix  is  sunpty  an 
averaged  effect  of  tte  ai^ular-dependent  reiadonriiip  be¬ 
tween  die  active  element  patterns  and  the  staini-alooe  ele^ 
ment  patterns.  Based  on  this  assumption,  various  ariiy 
calibration  inotxAnes  have  been  developed  in  vriikh  the 
ooiq£ng  matrix  is  detemimcd  by  mimmum  mean-sqoare  er- 
nn-  (MMSE)  mauhn^  of  the  two  sets  of  patterns  at  a  few 
known  incident  angles  [8;  9X 

While  both  of  these  apfnroadies  have  been  used  with 
success  on  arr^  oonsistliig  of  simple  di^e  elements,  the 
assumptions  u^  in  their  derivation  and  the  limitattons  of 
each  approadi  have  not  been  examined  in  detafl  in  the 
literature.  In  dns  p^icr,  we  set  out  to  compare  the  two 
different  appmachgs  frar  dgff»rmmwg  coupli^  matrix.  We 
first  review  the  two  approaches,  and  dfeenss  the  l^itarions  of 
each  aqiproach  in  Section  IL  In  Section  Ul,  numerical  results 
are  generated  using  die  standard  numerkal  elcctrcmiagnetic 
solver  NEC  (9]  to  verify  cmbt  ohservations.  In  Section  IV,  we 
show  that  the  ooo|ding-niatrix  model  can  be  further  extenefed 
fior  more  complex  antenna  arrays,  and  an  escanifde  is  provided 
to  fflnstrate  the  effectiveness  of  the  extended  modeL 

IL  IIUTUAL-COUPUI^  MODOS 

Assrniw  that  an  arr^  with  N  elements  is  located  in  the 
jy-fdane,  and  reodves  a  plaiK  wave  incident  from  angle  ^  in 
die  azhnutfa  plane.  Let  us  consider  a  matrix  reladonshfy 
bdween  the  active  eleinem  patterns  and  the  stand-akme 
element  patterns  as 

(1) 

where  and  A,j^  are  the  active  and  stand-alone  element 
patterns,  re^peedvefy,  and  M  is  the  mutual  coupin^  matrir 
Affaeo  has  N  rows,  widi  each  row  containmg  the  stand-akme 
efanewt  patrem  feu'a  partfenlar  element  m  the  array.  For  the 
ith  row;  it  can  be  ie|Hesented  as 

1^)  (2) 

wterc  K  the  stand-alone  etement  response  as  a  fime- 
thm  of  and  die  cxponemial  term  is  the  eaqplidt  position- 
dependent  phase  delay  ftir  the  /di  dement  located  at  (x^y^X 
The  stand-alone  pentern  matrix  of  the  array  is  then  as 

(3) ' 

The  active  {XRtem  matrix  has  the  strnctuie  as 
Adbeo*  excqpt  that  is  replaced  by  die  active  element 
pattern  of  the  f  th  demenL  The  two  different  approaches  m 
determining  the  nratoal  oonplii^  matrix  M  arc  discussed 
below. 
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A.  The  Z  Approach.  In  the  first  approach,  the  antenna  array 
with  N  elements  is  viewed  as  an  A^-port  microwave  network 
[3],  The  relationship  between  the  antenna  terminal  voltages 
and  currents  on  reception  under  a  plane-wave  excitation  can 
be  written  as 

Kj  “  ^\\^\  ^l.oc 

]/■  =  2/j/j  "f"  ***  ***  ^i,oc 

Vf^j  =  Zisf\I\  +  ••*  +  ***  '^Zf^f^If^  +  ^,oc* 

In  the  above  expression,  is  the  received  voltage  at  the 

/th  port  when  all  of  the  antenna  terminals  are  open  circuited, 
and  is  the  mutual  impedance  between  the  ith  and  the  ;th 
port  defined  as  =  k;//y|y^=o.  /t  If  we  assume  that  each 
antenna  is  loaded  by  load  impedance  Zy^  at  its  terminals  and 
use  the  relationship  Vi  =  -Z^^Ii,  we  can  express  the  termi¬ 
nal  voltages  in  the  following  matrix  form: 


or 

ZV  ==  V,,.  (5) 

Note  that,  if  we  stack  columnwise  the  received  voltage  vec¬ 
tors  V  obtained  for  different  plane-wave  incident  angles,  the 
resulting  matrix  is,  in  fact,  the  active  pattern  matrix  in 
Eq.  (1)  (apart  from  an  unimportant  scaling  factor).  Similarly, 
if  we  stack  the  open  circuit  voltage  vectors  and  assume 
that  they  are  the  same  as  the  voltages  received  under  stand¬ 
alone  conditions,  the  resulting  matrix  is  the  stand-alone  pat¬ 
tern  matrix  Therefore,  the  mutual  coupling  matrix  in 

(1)  can  be  written  as 

M  =  Z~'.  (6) 

This  completes  the  derivation  of  M  under  the  Z  approach. 

Since  the  normalized  impedance  Z  is  intrinsic  to  the 
structure  and  is  completely  independent  of  angle,  M  is  also 
angle  independent.  This  is  a  highly  desirable  conclusion. 
However,  it  hinges  on  the  assumption  that  the  received 
voltage  at  the  terminals  of  an  antenna  when  all  of  the 
antenna  elements  are  open  circuited  is  the  same  as  the 
received  voltage  when  all  of  the  other  elements  are  absent. 
From  the  electromagnetic  point  of  view,  this  is  clearly  not 


true  as  the  induced  currents  on  the  antenna  elements  will  not 
be  identically  zero  under  the  open-circuit  condition.  As  a 
result,  the  total  voltage  developed  at  the  antenna  terminals  is 
not  only  due  to  the  incident  field,  but  also  to  fields  produced 
by  the  induced  currents  on  other  elements.  Therefore,  while 
Eq.  (5)  is  always  true,  the  coupling  matrix  M  given  by  (6)  does 
not  rigorously  relate  the  active  and  stand-alone  element 
patterns.  Under  certain  situations,  however,  this  model  may 
work  approximately.  For  example,  if  we  consider  an  array 
comprised  of  half-wave  dipole  elements,  the  induced  currents 
on  every  element  should  be  very  small  under  the  open-circuit 
condition,  and  the  Z  approach  may  be  quite  adequate.  This 
will  be  examined  numerically  in  Section  III. 

B.  The  C  Approach.  In  the  second  approach,  the  mutual-cou¬ 
pling  matrix  is  viewed  as  an  average  of  the  angular-dependent 
relationship  between  the  active  and  stand-alone  element  pat¬ 
terns  [7].  Under  this  assumption,  M  can  be  obtained  by  using 
the  MMSE  matching  between  the  two  patterns  at  a  number 
of  incident  angles: 

\{  C!  ==  .  (7) 

We  refer  to  the  above  formulation  as  the  C  approach. 

The  effectiveness  of  the  C  approach  will  now  be  examined 
from  the  induced  current  point  of  view  based  on  the  method 
of  moments  (MoM).  In  MoM,  the  integral  equation  of  the 
induced  currents  on  ail  of  the  antennas  in  the  array  can  be 
written  in  a  matrix  form  as  follows: 

U  =  V  (8) 

where  I  is  the  current,  V  is  the  incident  plane-wave  excita¬ 
tion,  and  L  is  the  moment  matrix.  Suppose  that  we  remove  all 
but  the  /th  antenna  in  the  array;  we  can  find  the  stand-alone 
current  distribution  on  the  /th  antenna  by  solving  a  smaller 
moment  system  described  by 

-  y  .  (9) 

N  such  systems  can  be  written  for  all  of  the  elements  in  the 
array.  Since  the  right-hand  side  is  determined  by  the  incident 
field  only,  it  does  not  change  whether  the  antenna  is  stand¬ 
alone  or  in  the  array  environment.  Thus,  we  rewrite  Eq.  (8)  as 


I,  “ 

o 

1 _ 

I? 

IL] 

I 

= 

* 

>- 

0  L'J^ 

or 

LI  =  (10) 

By  solving  (10),  the  actual  current  distribution  is  related  to 
the  stand-alone  current  by 

(11) 

Next,  we  try  to  extract  the  currents  at  the  terminals  of  the 
antenna  elements  from  the  I  vector,  and  stack  them  over  all 
of  the  incident  angles  to  arrive  at  the  active  pattern  matrix 
A, rue*  Similarly,  we  try  to  extract  the  terminal  currents  from 
the  vector,  and  stack  them  to  form  the  stand-alone  pattern 
matrix  First,  we  sift  out  the  terminal  currents  from  I 
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and  the  corresponding  rows  from  the  matrix  (L~*L*^)  in  (11) 
to  form 

Itern.  =  (L''L")ltermI°  (12) 


where  the  subscript  “term”  indicates  the  rows  corresponding 
to  the  antenna  terminals.  Next,  we  sift  out  the  terminal 
currents  from  by  writing 


S, 


I«  = 


0 


0 


I 


0 

l.tcrm 


^W.ierm 


(13) 


where  S,  is  the  shape  of  the  current  distribution  on  the  /th 
element  and  is  its  corresponding  terminal  current. 

Substituting  (13)  into  (12),  we  can  find  the  relationship  be¬ 
tween  the  actual  and  stand-alone  terminal  currents.  Note 
that  the  moment  matrices  L  and  L®  are  angular  independent. 
If  the  shape  of  the  current  distribution  on  a  stand-alone 
antenna  element  is  also  angular  independent,  then  Eq.  (12) 
can  be  rigorously  cast  into  the  form  of  (1).  To  summarize,  we 
have  shown  that  the  relationship  between  the  active  element 
pattern  and  the  stand-alone  element  pattern  is  angular  inde¬ 
pendent  provided  that  the  shape  of  the  current  distribution 
on  each  stand-alone  antenna  element  is  independent  of  an¬ 
gle.  When  this  condition  is  satisfied,  the  mutual-coupling 
matrix  can  be  effectively  determined  by  (7). 

Let  us  now  examine  when  the  angular-independent  cur¬ 
rent  shape  condition  can  be  met.  The  first  situation  has  been 
discussed  in  [5],  where  all  of  the  antennas  are  vertical  wires 
and  the  incident  directions  have  the  same  elevation  angle.  In 
this  case,  the  stand-alone  current  distribution  of  an  antenna 
is  exactly  the  same  over  all  of  the  incident  angles.  A  second 
common  situation  is  when  the  antenna  elements  are  working 
near  resonance.  In  this  case,  the  stand-alone  current  distribu¬ 
tion  on  an  antenna  is  dominated  by  the  resonant  mode,  and 
thus  has  approximately  the  same  shape  over  all  of  the  inci¬ 
dent  angles.  An  example  of  this  situation  is  an  array  of 
half-wave  dipoles. 

III.  NUMERICAL  EXAMPLES 

In  this  section,  we  illustrate  the  conditions  for  the  Z  and  C 
approaches  to  be  valid  using  computer  simulation.  We  study 
a  linear  array  consists  of  four  dipoles,  whose  centers  are 
spaced  0.45  wavelength  apart,  as  shown  in  Figure  1.  The 
simulations  are  carried  out  using  NEC. 

In  the  first  example,  the  dipoles  are  a  half  wavelength  in 
length.  All  of  the  dipoles  are  loaded  with  73  ft  of  load 
impedance  at  the  center.  Elements  1  and  4  are  rotated  45 
and  —45®,  respectively,  about  the  y-axis.  Elements  2  and  3 
are  rotated  45  and  —45®,  respectively,  about  the  x-axis.  We 
determine  the  coupling  matrix  using  both  the  C  approach  and 
the  Z  approach  for  this  structure.  The  current  at  the  center 
of  each  antenna  is  computed  for  incident  angles  from  —  90  to 
90®  at  a  step  of  1®.  This  serves  as  the  reference  active  element 
pattern.  The  values  at  ±  75,  ±45,  and  0®  are  used  to  calculate 
the  C  matrix  from  (7).  To  determine  the  Z  matrix,  we  first 
obtain  the  mutual  admittance  by  short  circuiting  all  of  the 
antenna  ports,  driving  the  yth  antenna  with  a  voltage  source, 
and  computing  r  inverting  the  result¬ 

ing  generalized  admittance  matrix,  we  get  the  mutual 
impedance  Z^y.  The  normalized  impedance  matrix  Z  is  then 
obtained  from  (5). 


Once  we  obtain  the  two  coupling  matrices,  we  use  them  to 
calculate  the  active  element  patterns  from  the  stand-alone 
element  patterns.  Notice  that,  since  the  elements  are  tilted, 
even  the  stand-alone  element  patterns  are  not  isotropic.  The 
active  element  patterns  of  antenna  elements  1  and  2  are 
plotted  in  Figure  2(a)  and  (b),  respectively.  The  solid  curves 
are  the  active  element  patterns  directly  computed  by  NEC. 
The  dashed  curves  are  obtained  using  the  Z  approach,  and 
the  dotted  lines  are  modeled  using  the  C  approach.  We  find 
that  the  dashed  curves  are  quite  close  to  the  computed  ones, 
but  are  not  exactly  the  same.  This  is  because  the  condi¬ 
tion  described  in  Section  II-A  is  only  approximately  satisfied, 
even  for  an  array  of  half-wave  dipoles.  On  the  other  hand, 
the  curves  using  the  C  approach  match  almost  perfectly  with 
the  computed  ones.  In  this  case,  the  antennas  operate  in 
resonance,  and  the  stand-alone  current  distribution  has  the 
same  shape  over  all  incident  angles,  even  though  they  are  not 
all  vertically  aligned. 

To  further  observe  the  effectiveness  of  the  mutual-cou¬ 
pling  models,  we  apply  the  active  element  patterns  to  a 
direction-finding  problem,  which  is  quite  sensitive  to  the 
element  patterns.  Two  uncorrelated  incoming  signals  are 
assumed  from  the  incident  angles  of  50  and  65®.  The  signal- 
to-noise  ratio  (SNR)  is  set  to  30  dB.  Tlie  angles  of  arrival  are 
determined  using  the  standard  MUSIC  algorithm,  and  the 
resulting  normalized  power  spectra  are  plotted  in  Figure  3. 
The  solid  curve  is  the  result  from  the  reference  active  ele¬ 
ment  patterns.  For  comparison,  the  result  by  neglecting  the 
mutual-coupling  effect  is  obtained  by  using  the  stand-alone 
element  patterns,  and  is  plotted  as  the  dash-dot  curve.  The 
large  degradation  shows  the  importance  of  taking  mutual 
coupling  into  consideration.  The  dashed  curve  is  the  result  of 
the  Z  approach.  The  dotted  curve  is  obtained  using  the  C 
approach.  It  is  observed  that  the  C  approach  performs  almost 
as  well  as  the  reference  active  element  patterns,  while  the  Z 
approach  leads  to  some  degradation  in  the  peak  position  and 
width  of  the  power  spectra. 

In  the  second  example,  we  replace  the  half-wave  dipole 
elements  with  full-wave  dipoles.  Compared  with  the  previous 
example,  the  condition  is  strongly  violated  for  the  Z 
approach,  and  the  near-resonant  condition  is  not  satisfied  for 
the  C  approach.  Figure  4  shows  the  computed  and  modeled 
active  element  patterns.  This  time,  neither  the  Z  nor  the  C 
approach  is  close  to  the  reference  result. 

Finally,  we  make  the  full-wave  dipoles  perpendicular  to 
the  jiy-plane.^Now,  the  vertical  wire  condition  is  satisfied  for 
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(a)  clement  1 


(b)  element  2 

Figure  2  Computed  and  modeled  active  element  patterns  for  array  of  rotated  half-wave  dipoles 
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the  C  approach.  The  resulting  element  patterns  are  plotted  in 
Figure  5.  The  Z  approach  still  works  poorly  as  in  example  2. 
However,  the  C  approach  perfectly  models  the  active  element 
patterns.  As  discussed  in  Section  II-B,  the  matching  in  this 
case  is  exact. 

IV  OP  jHE  mutual-coupling  MODEL 

As  can  be  seen  in  the  examples,  the  C  approach  is  more 
accurate  than  the  Z  approach.  Yet,  the  situations  when  the  C 
approach  can  be  applied  are  still  quite  limited.  The  approach 
fails  for  more  complex  antenna  structures.  By  observing  Eq. 
(11),  we  can  see  that,  if  the  array  elements  consist  of  thin 
wires,  each  piece  of  wire  affects  the  active  element  pattern  in 
the  same  manner  as  a  dipole  antenna.  Thus,  for  certain  cases, 
the  model  in  (1)  can  be  extended  to  include  the  coupling 
from  all  parts  of  the  array.  For  example,  if  we  have  an  array 
of  Yagi  antennas,  we  may  include  the  coupling  from  the 
parasitic  elements  to  model  the  active  element  patterns  more 
accurately.  More  specifically,  we  extend  the  model  as  follows: 


A 


true 


[c 


C'] 


^theo 

^iheo 


(14) 


where  A\hco  the  stand-alone  receiving  patterns  of  the 
parasitic  elements  and  C'  describes  the  coupling  from  the 


(b)  Yagi  element 


90 


Figure  6  Four-element  array  with  Yagi  elements 


parasitic  elements  to  the  driven  elements.  Like  the  standard 
C  approach,  (14)  is  valid  if  all  of  the  parasitic  elements  are 
vertical  or  their  lengths  are  close  to  the  resonant  length.  The 
MMSE  matching  can  again  be  used  to  determine  the  cou¬ 
pling  matrix,  except  that  the  number  of  incident  angles  needed 
to  give  a  unique  solution  is  increased. 

As  an  example,  we  look  at  an  array  consisting  of  four 
parallel  Yagi  antennas,  as  shown  in  Figure  6(a).  All  of  the 
elements  in  the  array  are  perpendicular  to  the  jy-plane,  and 
the  separation  between  two  adjacent  Yagis  is  one  wavelength. 
Each  Yagi  antenna  consists  of  three  wires,  with  the  center 
wire  being  the  driven  element,  as  shown  in  Figure  6(b).  The 
stand-alone  element  pattern  is  plotted  in  Figure  6(c).  Next, 
the  coupling  matrix  in  (14)  is  determined  from  the  computed 


90 


(b)  element  2 


Figure  7  Computed  and  modeled  active  element  patterns  for  an 
array  of  Yagi  antennas 
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active  element  patterns  at  12  incident  angles  uniformly  dis¬ 
tributed  between  —180  and  180®.  The  reference  and  the 
modeled  active  element  patterns  of  the  first  and  second 
antennas  are  plotted  in  Figure  7(a)  and  (b),  respectively.  The 
solid  curves  are  the  reference  patterns^  and  the  dotted  curves 
are  obtained  using  the  extended-C  approach.  They  match 
with  the  reference  patterns  very  well.  As  a  comparison,  we 
also  model  the  active  element  patterns  using  the  Z  approach 
and  the  C  approach.  The  results  are  shown  in  dashes  and 
dash“dot  curves,  respectively.  The  agreement  is  obviously  not 
as  good.  This  is  because  the  open-circuit  voltage  condition 
does  not  hold  for  the  Z  approach,  and  the  constant  current 
shape  condition  for  each  array  element  is  violated  for  the  C 
approach. 

V*  CONCLUSIONS 

In  this  paper,  the  two  commonly  used  approaches  for  deter¬ 
mining  the  mutual-coupling  matrix  in  array  antennas  have 
been  discussed.  It  has  t^en  shown  that  the  Z  approach  is  an 
approximation  that  applies  when  the  element  pattern  under 
the  open-circuit  condition  can  be  well  approximated  by  the 
stand-alone  element  pattern.  It  has  also  been  shown  that  the 
C  approach  holds  true  whenever  the  current  distribution  of  a 
stand-alone  element  has  the  same  shape  over  all  incident 
angles.  This  condition  is  satisfied  for  vertical  wire  antennas  or 
near-resonant  antennas.  We  have  found  from  our  numerical 
examples  that  the  C  approach  is,  in  general,  more  accurate 
than  the  Z  approach.  We  have  also  extended  the  C  approach 
to  model  more  complex  wire  elements  by  considering  each 
wire  as  an  individual  antenna.  An  example  of  an  Yagi  array 
has  been  provided  to  demonstrate  the  effectiveness  of  the 
extended-C  model. 
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Abstract — A  frequency  extrapolation  technique  is  proposed  to  obtain 
the  broad  band  radiation  patterns  of  antenna-platform  radiation 
problems.  A  frequency-dependent  time-of-arrival  model  is  applied 
to  the  induced  current  computed  at  low  frequencies.  The  model 
parameters  are  estimated  by  applying  the  ESPRIT  superresolution 
algorithm  to  the  computed  data.  A  pre-multiplication  scheme  in 
conjunction  with  the  complex  time-of-  arrival  estimation  from  ESPRIT 
is  used  to  determine  the  additional  frequency-dependent  factors  in  the 
model.  The  current  at  high  frequencies  is  then  extrapolated  based  on 
the  model  and  the  reidiated  field  is  computed  using  the  extrapolated 
current.  The  algorithm  is  applied  to  several  2D  and  3D  antenna- 
platform  radiation  problems.  The  extrapolated  radiation  patterns  axe 
found  to  be  in  good  agreement  with  those  generated  by  brute  force 
computation.  Since  the  current  required  for  the  extrapolation  is  only 
computed  at  lower  frequencies,  large  savings  in  computation  time  and 
memory  can  be  achieved. 

1  Introduction 

2  EVequency  Dependent  Model  and  Determination  of 
Model  Parameters 

3  Algorithm  Testing  and  Error  Estimation  Using  Simu¬ 
lated  Data 

4  Numerical  Results 

5  Conclusions 
Acknowledgment 
References 


866 


Su,  Wang,  and  Ling 


1.  INTRODUCTION 

The  numerical  solution  of  electromagnetic  scattering  and  radiation 
problems  from  complex  structures  is  usually  very  computer  intensive 
in  time  and  memory,  especially  at  high  frequencies.  In  the  standard 
method  of  moments  (MoM),  for  example,  the  computation  complexity 
scales  as  where  /  is  the  frequency.  Even  for  fast  solvers  such  as 
the  fast  multipole  method  [1,  2],  the  frequency  scaling  is  still  between 
P  and  /^.  Thus  it  is  desirable  to  obtain  the  broad  band  response 
from  the  computation  results  at  only  a  few  frequencies.  This  can  be 
accomplished  by  fitting  the  computed  frequency  response  to  a  reduced- 
order  model  and  extracting  the  model  parameters  from  the  data. 
Approaches  based  on  model-based  parameter  estimation  have  been 
studied  extensively  and  applied  to  many  aspects  in  electromagnetic 
computation  [3-11].  The  appropriate  reduced-order  models  to 
parameterize  the  physical  observables  are  different  depending  on  the 
frequency  region  of  interest.  In  the  low  frequency  range,  the  rational 
function  model  is  widely  used,  while  in  the  high  frequency  range,  the 
time-of-arrival  model  is  the  preferred  choice  in  characterizing  the  ray- 
optical  behaviors  of  fields  and  currents. 

In  this  paper,  we  address  the  model-based  frequency  extrapolation 
of  antenna-platform  radiation  problems.  It  is  well  known  that  the 
radiation  patterns  of  antennas  are  strongly  affected  by  the  mounting 
platform.  Since  the  simulation  of  such  radiation  problems  is  very 
costly  to  generate  as  a  function  of  frequency,  we  set  out  to  extrapolate 
the  high  frequency  response  from  computed  data  at  lower  frequencies. 
Our  approach  entails  fitting  the  computed  current  at  low  frequencies 
to  a  time-of-arrival  model  and  estimating  the  model  parameters  by 
the  superresolution  algorithm  ESPRIT  [12].  Previously,  this  approach 
has  been  applied  to  radar  signature  prediction  with  good  success 
[10].  In  that  work,  the  coefficients  of  the  time-of-arrival  model  were 
assumed  to  be  frequency  independent.  In  the  radiation  problem,  the 
primary  radiation  from  the  antenna  is  usually  frequency  dependent. 
Furthermore,  when  there  exist  higher-order  multiple  interactions  on 
the  platform,  the  amplitude  of  each  mechanism  is  in  general  also 
frequency  dependent  [13-17].  Thus  a  more  accurate  model  is  needed 
to  parameterize  the  current.  Here  we  generalize  the  extrapolation 
algorithm  by  using  a  frequency-dependent  time-of-arrival  model.  To 
extract  the  additional  frequency  dependence  in  this  model,  we  adopt 
an  approach  similar  to  the  one  proposed  in  [15],  which  was  developed 
to  effectively  extract  the  frequency  dependency  of  the  scattering 
mechanisms  in  measured  backscattered  data.  Our  results  show  that 
the  radiation  pattern  can  be  accurately  extrapolated  based  on  the 
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frequency  dependent  model. 

This  paper  is  organized  as  follows.  Section  2  describes  the 
frequency-dependent  time-of-arrival  model  and  the  procedure  to 
determine  the  model  parameters.  Section  3  gives  a  numerical  example 
of  the  extrapolation  algorithm  using  simulated  data.  The  performance 
of  the  algorithm  in  the  presence  of  noise  is  investigated  and  the  errors 
in  the  estimation  of  model  parameters  are  quantified.  In  Section  4  we 
apply  the  algorithm  to  extrapolate  the  induced  current  in  antenna- 
platform  radiation  problems.  The  radiation  patterns  predicted  at  high 
frequencies  are  compared  against  the  reference  MoM  computation  in 
both  2D  and  3D  platforms.  Conclusions  are  given  in  Section  5. 


Figure  1.  Time-of-arrival  model  for  the  induced  current  at  point  P 
accounts  for  the  direct  incident  radiation  from  the  antenna  and  the 
multiple  scattered  waves  from  other  parts  of  the  platform. 


2.  FREQUENCY  DEPENDENT  MODEL  AND 
DETERMINATION  OF  MODEL  PARAMETERS 

We  postulate  that  the  current  induced  on  a  complex  platform  due  to 
illumination  from  a  primary  source  can  be  well  described  by  a  time-of- 
arrival  model  at  high  frequencies.  The  different  time-of-arrival  terms 
correspond  to  the  different  incident  wave  mechanisms  from  both  the 
direct  antenna  radiation  and  higher-order  scattering  from  other  parts 
of  the  platform,  as  illustrated  in  Fig.  1.  Therefore,  the  current  at  an 
arbitrary  point  P  on  the  platform  can  be  expressed  as: 

=  (1) 

n 

where  tn  is  the  arrival  time  of  the  nth  incident  wave  and  an  is  its 
amplitude.  In  [10],  a„  was  assumed  to  be  independent  of  frequency. 
However,  in  the  antenna-platform  radiation  problem,  the  primary 
radiation  from  the  antenna  is  in  general  frequency  dependent.  For 


868 


Su,  Wang,  and  Ling 


instance,  the  radiation  strength  of  a  dipole  antenna  is  proportional  to 
the  u)  (or  for  a  line  source  in  2D  problems).  Furthermore,  for 
complex  platforms,  the  incident  waves  to  a  specific  current  element 
can  come  from  the  scattered  fields  from  other  parts  of  the  platform. 
These  secondary  incident  fields  are  also  frequency  dependent.  For 
canonical  shapes,  their  exact  frequency  dependencies  are  predicted  by 
the  geometrical  theory  of  diffraction  (GTD),  and  are  in  the  form  of 
(jj^  where  ex.  is  the  frequency  factor.  For  example,  the  strength  of 
the  scattered  field  from  corner  diffraction  has  a  frequency  dependency 
of  and  that  from  edge  diffraction  has  a  frequency  dependency 
of  Therefore,  the  time-of-arrival  model  can  be  improved  by 

incorporating  the  frequency  factor: 

=  (2) 

n 

where  an,  is  the  frequency  factor  for  each  incident  wave.  Compared 
with  (1),  this  model  is  better  matched  to  the  actual  physics  of  the 
radiation  and  scattering  mechanisms. 

Since  we  are  attempting  to  extrapolate  the  frequency  response 
to  a  much  broader  range,  the  accurate  estimation  of  the  frequency 
faotors  is  critical.  A  small  error  in  a  will  result  in  dramatic  difference 
in  amplitude  at  frequencies  in  the  extrapolated  region.  However,  the 
existing  superresolution  algorithms  based  on  eigenspace  decomposition 
(e.g.,  ESPRIT  and  MUSIC)  apply  only  to  signals  in  the  form  of  (1).  For 
instance,  the  ESPRIT  algorithm  [12]  estimates  (an,^n)  from  uniformly 
sampled  current  data  from  cji,  to  ujm  based  on  the  data  model: 

^  =  ^<^1,^2, . . . ,  wm  (3) 

n 

where  n  denotes  additive  Gaussian  white  noise.  This  model  can 
be  easily  extended  to  estimate  complex  time-of-arrivals  via  analytic 
continuation.  If  we  separate  the  real  and  imaginary  part  of  tn  the 
resulting  model  can  be  written  as 

J{u)  =  Y.  (4) 

n 

where  the  tilda  symbol  is  used  to  indicate  a  complex  number.  This 
is  the  well-known  complex  exponential  model  and  has  been  used  to 
approximately  model  the  frequency  dependence  in  a„  [13].  Comparing 
(4)  to  the  frequency-dependent  time-of-arrival  model  (2),  we  find  that 
the  only  difference  between  the  two  is  the  frequency  dependency,  which 
is  in  the  form  of  a;“  in  (2)  and  in  (4).  For  a  narrow  band  signal, 
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(4)  is  usually  a  good  approximation  of  (2),  and  can  be  simply  used  to 
model  the  frequency  dependence  within  this  band.  However,  in  order  to 
achieve  accurate  extrapolation  results,  the  complex  exponential  model 
is  a  poor  choice  since  it  is  not  properly  matched  to  the  underlying 
physics. 

To  better  determine  the  frequency  factor  an,  we  adopt  an 
approach  similar  to  the  one  proposed  in  [15]  based  on  a  pre¬ 
multiplication  concept.  We  multiply  the  data  by  an  assumed  frequency 
factor  a;-"'".  The  data  model  then  becomes 

J(w)  =  (5) 

n 

Next,  we  apply  ESPRIT  and  the  complex  exponential  model  on  J{u)) 
to  estimate  the  real  and  imaginary  parts  of  in-  Note  that  if  the 
frequency  dependence  of  the  pre-multiplier  is  exactly  matched  to  that 
of  a  particular  mechanism,  the  resulting  exponential  term  will  have  no 
frequency  dependence.  Consequently,  the  estimated  complex  time-of- 
arrival  term  will  have  a  zero  imaginary  part.  Therefore,  by  repeatedly 
pre-multiplying  the  signal  using  different  values  of  am  and  observing 
the  imaginary  parts  of  the  resulting  in,  we  can  detect  the  correct 
frequency  dependence  whenever  Im(f„)  goes  to  zero.  The  implicit 
assumption  of  this  approach  is  that  the  mismatched  terms  will  not 
significantly  distort  the  estimation  of  Im(f„)  of  the  matched  term. 
This  is  usually  true  when  the  arrival  times  Re(tn)  are  well  separated. 
When  two  or  more  of  the  arrival  times  are  very  close  to  each  other, 
this  algorithm  becomes  less  accurate  due  to  the  interference  between 
the  exponential  terms.  The  imaginary  part  of  may  not  vanish 
even  when  the  corresponding  frequency  factor  is  matched  by  the  pre¬ 
multiplication.  Quantitative  evaluation  of  the  resulting  error  will  be 
discussed  in  detail  in  the  next  section.  Once  an  and  tn  are  estimated, 
the  amplitude  an  can  be  easily  found  by  the  standard  least  squares 
procedure. 

3.  ALGORITHM  TESTING  AND  ERROR  ESTIMATION 
USING  SIMULATED  DATA 

To  test  the  extrapolation  algorithm,  we  consider  the  signal 

J{(jj)  =  (1  +  -I-  (4.3  —  +  n, 

a;  =  5,6, . . . ,35  (6) 

where  n  is  additive  Gaussian  white  noise.  The  signal-to-noise  ratio 
(SNR)  is  set  to  10  dB.  We  now  use  a  two-term  time-of-arrival  model 
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Figure  2.  Extrapolation  of  a  simulated  signal  using  different  models. 

to  extrapolate  the  signal  based  on  its  first  10  samples.  The  actual 
signal  is  plotted  as  the  solid  line  in  Fig.  2.  The  extrapolation  results 
of  three  different  models  are  plotted.  The  frequency  independent 
model,  which  is  plotted  in  dashed  dot,  is  obtained  by  the  model  in 
(1)  with  real  time-of-ar rival  tn  determined  by  ESPRIT.  This  model 
matches  badly  with  the  original  curve,  even  in  the  sampled  region.  The 
complex  exponential  model,  plotted  in  dashed  line,  uses  the  model  in 
(4)  with  complex  in-  The  resulting  curve  is  in  good  agreement  with  the 
actual  signal  within  the  sampled  region,  but  not  in  the  extrapolated 
region.  This  indicates  that  although  the  model  (4)  can  be  a  very  good 
approximation  to  the  actual  signal  in  the  sampled  region,  it  cannot 
be  used  to  accurately  extrapolate  the  signal  because  it  does  not  have 
the  correct  frequency  dependency.  The  frequency-dependent  time-of- 
arrival  model,  which  is  plotted  in  dots,  matches  the  best  with  the 
original  signal  in  both  the  sampled  region  and  the  extrapolated  region. 
Its  parameters  are  estimated  using  the  method  described  in  Section  2. 

To  show  the  detail  workings  of  the  pre-multiplier  procedure,  the 
imaginary  parts  of  tn  corresponding  to  the  two  terms  in  (6)  are  plotted 
as  functions  of  Om  in  Figs.  3(a)  and  3(b)  at  the  SNR  levels  of  oo  and 
10  dB,  respectively.  The  two  curves  cross  zero  at  the  corresponding 
frequency  factor  of  —0.5  and  0.5.  It  can  be  observed  that  the  shapes 
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Figure  3.  Variation  of  Im(tn)  as  a  function  of  the  pre-multiplication 
factor  am  for  the  simulated  signal  at  two  different  SNR  levels. 


of  the  curves  change  when  the  SNR  is  decreased,  especially  for  the 
term  with  smaller  power  (a„  =  -0.5).  However,  the  positions  of  the 
zero-crossing  points  vary  only  weakly  as  the  SNR  is  decreased. 

Next,  we  examine  the  error  performance  of  the  extrapolation 
algorithm  through  numerical  simulation.  Among  the  three  sets  of 
model  parameters  a„,  a„  and  the  time-of-arrival  is  directly 
estimated  by  the  ESPRIT  algorithm.  When  frequency  dependent 
mechanisms  are  present  in  the  signal,  the  estimation  error  in  tn  is 
expected  to  be  larger  due  to  the  interference  between  the  exponential 
terms,  especially  when  two  or  more  arrival  times  are  very  close  to  each 
other.  In  this  example,  we  estimate  the  arrival  times  in  a  simulated 
signal  consisting  of  two  frequency  dependent  exponential  terms  using 
ESPRIT.  As  a  comparison,  we  also  determine  the  arrival  times  in  a 
two-term  frequency  independent  signal.  In  both  signals,  the  power 
of  the  stronger  term  is  five  times  that  of  the  weaker  term.  For  the 
frequency  dependent  signal,  the  two  frequency  factors  are  set  to  0.5 
and  —0.5  for  the  stronger  and  weaker  terms,  respectively.  The  test 
is  performed  1000  times  subject  to  random  noise.  In  Fig.  4(a),  the 
root  mean  squared  (RMS)  error  on  the  time-of-arrival  estimate  of  the 
weaker  term  is  plotted  as  a  function  of  the  separation  between  the  two 
arrival  times  and  for  different  SNR  conditions.  The  dashed  lines  are 
the  estimation  errors  from  the  frequency  independent  signal,  while  the 
solid  lines  are  the  errors  from  the  frequency  dependent  signal.  Both 
axes  are  calibrated  in  terms  of  the  Fourier  resolution,  which  is  the 
reciprocal  of  the  available  data  bandwidth.  It  is  shown  that  the  arrival 
time  error  is  larger  for  the  frequency  dependent  signal,  especially  when 
the  SNR  is  high.  However,  the  degradation  is  not  too  severe  overall. 


872  Su,  Wang,  and  Ling 

Good  estimation  (with  error  less  than  0.1  Fourier  bin)  is  still  achievable 
for  signal  separation  well  within  one  Fourier  bin. 

The  RMS  error  of  estimating  the  frequency  factor  an  of  the  weaker 
signal  is  plotted  in  Fig.  4(b)  for  the  frequency  dependent  signals  used  in 
the  above  test.  We  observe  that  the  error  decreases  as  the  separation 
between  the  exponential  terms  increases  and  as  the  SNR  increases, 
similar  to  the  trends  in  Fig.  4(a).  When  the  two  arrival  times  are 
close,  the  interaction  between  the  two  terms  results  in  error  on  the 
estimation  of  and  the  error  is  strongly  affected  by  the  noise  level. 
We  can  further  decrease  the  error  by  imposing  the  constraint  that  the 
frequency  factors  are  integer  multiples  of  1/2,  as  dictated  by  GTD. 
This  is  implemented  when  we  deal  with  actual  electromagnetic  data  in 
the  next  section. 

Finally  we  look  at  the  extrapolation  error  of  the  entire  frequency- 
dependent  time-of-arrival  model.  An  extrapolation  ratio  of  3  to  1  is 
used  in  the  example,  i.e.,  the  frequency  response  is  extrapolated  to 
three  times  the  original  bandwidth.  The  percentage  RMS  error  is 
plotted  in  Fig.  4(c).  When  the  separation  between  the  two  terms 
is  large,  the  extrapolation  error  behaves  similar  to  the  errors  of  tn 
and  an-  When  the  two  terms  become  too  close  in  time,  we  note  that 
the  error  curves  actually  reach  a  maximum  and  then  drop  off.  This 
is  because  below  this  point,  the  two  exponential  terms  are  too  close 
to  be  resolved  within  the  extrapolation  region  of  interest.  For  larger 
extrapolation  bandwidth,  the  position  of  the  peak  should  move  toward 
zero. 

4.  NUMERICAL  RESULTS 

We  now  apply  the  frequency  extrapolation  technique  to  computation 
data  from  antenna-platform  radiation  problems.  In  the  first  example,  a 
two-dimensional  structure  shown  in  Fig.  5  is  considered.  The  platform 
is  14  m  in  length  and  3  m  in  height.  The  antenna  is  a  horizontal 
line  source  placed  at  5  m  above  the  platform.  The  induced  current 
on  the  platform  is  computed  from  0.1  to  0.5  GHz  at  21  frequency 
points  using  2D  MoM.  The  current  is  extrapolated  to  1.3  GHz  and 
radiated  field  is  then  computed  based  on  the  extrapolated  current. 
Both  the  frequency-independent  and  frequency-dependent  time-of- 
arrival  models  axe  used  to  perform  the  extrapolation,  and  the  resulting 
radiated  fields  at  the  observation  angle  oiO  =  40°  are  plotted  in  Figs. 
6(a)  and  6(b),  respectively.  Also  plotted  is  the  reference  MoM  results 
obtained  via  brute  force  computation.  The  primary  radiation  of  the 
dipole  antenna  is  not  included  in  the  plots  so  that  we  can  better 
observe  the  secondary  radiation  from  the  platform.  It  is  obvious 
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(a)  RMS  error  in  estimating  the  time-of-arrival  tn  using  ESPRIT. 


(b)  RMS  error  in  estimating  the  frequency  factor 
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Separation  between  two  arrival  times  (1/BW) 

(c)  Percentage  RMS  error  of  the  extrapolated  signal  at 
the  extrapolation  ratio  of  3:1. 

Figure  4.  Error  performance  of  the  extrapolation  algorithm  as  a 
function  of  the  separation  between  two  arrival  times  at  different  SNR 
levels. 


I 
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Figure  5.  2D  platform  geometry. 

that  the  frequency  dependency  in  the  field  response  is  not  captured 
by  the  frequency-independent  time-of-arrival  model,  while  the  field 
predicted  by  the  frequency  dependent  model  is  in  good  agreement 
with  the  computed  result.  The  time  domain  response  corresponding 
to  Fig.  6(b)  is  plotted  in  Fig.  6(c).  It  is  shown  that  the  time-domain 
peaks  are  well  characterized  by  the  extrapolated  field.  Finally,  the 
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(c)  Time  domain  response  predicted  by  the  frequency 
dependent  model. 

Figure  6.  Frequency  extrapolation  example  for  the  2D  platform. 


radiation  patterns  obtained  through  brute  force  MoM  computation 
and  frequency  extrapolation  are  plotted  as  functions  of  frequency  and 
observation  angles  in  Figs.  7(a)  and  7(b),  respectively.  Very  good 
agreement  is  observed  between  the  two  figures.  For  a  quantitative 
comparison,  the  correlation  index  between  the  two  figures  is  computed 
using  the  definition 


//  Etif,e)E2{f,e)dfde 

R=TTrr - ^ 

iyy  (\Ei{f,e)\''  +  \E2{fM)dfd9 

where  Ei  and  E2  are  the  radiated  fields  in  linear  scale  obtained  by 
extrapolation  and  MoM  computation,  respectively.  The  correlation 
index  between  the  two  figures  is  0.9992  in  the  sampled  region  and 
0.9857  in  the  extrapolated  region.  As  expected,  the  correlation  is  lower 
in  the  extrapolated  region.  However,  the  drop  off  is  relatively  small. 

Next,  we  look  at  a  3D  platform  shown  in  Fig.  8.  The  antenna  is 
a  horizontal  dipole  oriented  in  the  x  direction.  The  induced  current  is 
computed  from  0.1  to  0.36  GHz  at  13  frequencies  and  extrapolated  to 
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Figure  7.  Comparison  of  the  radiated  field  generated  from  brute 
force  MoM  computation  and  frequency  extrapolation  as  a  function  of 
frequency  and  observation  angle. 
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Figure  8,  3D  platform  geometry. 


0.7  GHz.  The  computation  is  carried  out  using  FISC  [18],  which  is  a 
3D  MoM  code  based  on  the  fast  multipole  method.  The  extrapolated 
frequency  and  time  domain  responses  at  the  observation  angle  <f)ei  = 
30'',  ^az  =  -60"^  are  plotted  in  Figs.  9(a)  and  9(b),  respectively. 
Also  plotted  for  comparison  are  the  reference  responses  computed  by 
FISC  via  brute  force.  The  major  radiation  features  are  captured  in 
both  domains  by  the  extrapolation.  In  the  time  domain  response, 
we  observe  that  the  amplitudes  of  some  of  the  weaker  peaks  are  not 
correctly  predicted  by  the  extrapolation,  although  their  arrival  times 
are  estimated  quite  accurately.  We  believe  this  is  due  to  the  estimation 
error  of  an  in  the  current  parameterization,  as  we  expect  larger  errors 
in  the  frequency  factors  for  the  weaker  time-of-arrival  terms.  Figs. 
10(a)  and  10(b)  show  the  reference  and  extrapolated  radiation  patterns 
as  functions  of  frequency  and  azimuth  angle  when  the  elevation  angle  is 
fixed  at  50°.  Good  qualitative  agreement  is  observed.  The  correlation 
index  between  the  two  figures  is  found  to  be  0.9980  in  the  sampled 
region  and  0.9742  in  the  extrapolated  region,  respectively.  This  result 
is  a  little  worse  than  the  2D  example  because  the  noise  level  of  the 
FISC-computed  current  is  higher  than  that  of  the  2D  MoM  code.  As 
was  shown  in  the  last  section,  the  extrapolation  performance  is  affected 
by  the  SNR  in  the  computed  data.  The  computation  time  of  the  brute 
force  reference  results  is  about  50  hours  on  a  Pentiumll  400  MHz  PC, 
while  the  total  computation  time  to  carry  out  the  EM  computation  in 
the  sampled  region  and  the  extrapolation  process  is  7  hours. 


OO 
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(a)  Frequency  response  predicted  by  the  frequency 
dependent  model  at  ^el  =  30°,  (t>az  =  —60°. 


(b)  Time  domain  response  predicted  by  the 
frequency  dependent  model. 


Figure  9.  Frequency  extrapolation  example  for  the  3D  platform. 
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Figure  10.  Comparison  of  the  radiated  field  generated  from  brute 
force  FISC  computation  and  frequency  extrapolation  as  a  function  of 
frequency  and  azimuth  angle  at  the  elevation  angle  of  <j>ei  =  50°. 
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In  this  paper,  the  frequency-dependent  time-of-arrival  model  has  been 
applied  to  extrapolate  the  induced  current  and  radiation  pattern  in 
antenna-platform  radiation  problems.  The  model  parameters  are 
estimated  by  applying  the  ESPRIT  superresolution  algorithm  to  the 
computed  data.  A  pre-multiplication  scheme  in  conjunction  with  the 
complex  time-of-arrival  estimation  from  ESPRIT  is  used  to  determine 
the  additional  frequency-dependent  factors.  The  performance  of  the 
algorithm  in  the  presence  of  noise  has  been  evaluated  based  on 
simulated  data  and  errors  in  the  estimation  of  model  parameters  have 
been  quantified.  Our  results  show  that  the  method  is  quite  robust. 
The  algorithm  has  been  applied  to  extrapolate  the  induced  currents 
and  radiation  patterns  in  both  2D  and  3D  antenna-platform  radiation 
problems.  The  radiation  patterns  computed  from  the  extrapolated 
currents  have  been  found  to  be  in  good  agreement  with  those  generated 
by  brute  force  computation. 

Although  the  determination  of  model  parameters  is  more 
complicated,  the  frequency  dependent  model  show  significant 
performance  improvement  over  the  frequency  independent  model.  This 
is  due  to  the  improved  modeling  of  the  scattering  physics.  Since 
the  current  required  for  the  extrapolation  is  only  computed  at  lower 
frequencies,  large  savings  in  computation  time  and  memory  can  be 
achieved. 
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Abstract — This  paper  analyzes  a  set  of  first  order  triangular  patch 
basis  functions  fisr  the  method  of  moments  (MoM)  isolation  of 
electromagnetic  scattering  problems,  and  compares  its  pafbrmanoe 
to  the  traditional  Rao-A^ton-Glismn  (RWG)  basis.  These  bngig 
functions  provide  a  &ster  convergence  and  a  smootlnH-  snrferg  current 
representation.  It  is  shown  that  the  computational  complexity 
involved  is  about  the  same  as  the  RWG  basis,  which  mpang  that. 
existing  computer  codes  can  be  ad^ted  to  use  them  with  TniTiimal 
modifications. 
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1.  INTRODUCTION 

Since  the  introduction  of  the  Rao-Wilton-Glisson  (RWG)  triangular 
patdi  basis  functions  [1],  th^  have  become  the  most  widely  used 
basis  functions  for  solving  electromagnetic  scattering  problems  with  the 
method  of  moments.  These  basis  functions  have  the  highly  HpRirahlA 
property  of  being  &ee  of  line  diarges  while  being  able  to  model  any 
arbitrary  shape  surface.  For  this  reason  they  have  been  used  by  many 
authors  [2—9].  Improvements  in  the  original  formulation  such  as  using 
curvilinear  triangular  patches  [10]  and  incorporating  edge  mTiditions 
[11]  have  also  been  reported. 

Some  authors  [12]  have  pointed  out  that  the  results  obtained  naing 
the  RWG  basis  are  sometimes  very  dependent  on  the  triangulation 
scheme  used.  FVirthermore,  it  is  wdl  recognized  that  the  resulting 
current  distribution  often  contains  discontinuities  that  hinHAr  its  use 
for  accurate  current  and  near  field  prediction.  Filtering  scheme  has 
been  proposed  to  alleviate  this  problem  [13],  but  leads  to  loss  of  spatial 
resolution.  These  aforementioned  problems  are  due  to  the  fact  that 
the  RWG  basis  functions  cannot  represent  an  arbitrary  KnAar  current 
distribution  on  each  triangle. 

Ih  [14],  Wang  and  Wfebb  presented  a  new  set  of  hierarchical  basis 
functions  suitable  for  a  p-adaptation  sdieme  firom  linAar  to  bigW 
order.  They  show  that  much  better  results  can  be  obtained,  for  the 
same  meshing,  with  second  or  higher  order  basis  functions.  To  use 
these  bases,  however,  a  new  code  has  to  be  written. 

In  this  paper  we  present  an  enhanced  set  of  triangular  patch 
basis  functions,  which  retain  the  same  nice  properties  as  the  standard 
RWG  basis,  but  are  capable  of  represent  any  linear  current  distribution 
over  each  trianj^e.  With  the  improvement,  these  basis  fiitirtmTig  are 
much  less  sensitive  to  the  triangulation  sdieme  and  provide  a  mu**h 
better  representation  of  the  real  current  distribution,  leading  to  a 
faster  convergence  of  the  moment  method  solution.  These  fiiTirtin^ip  are 
equivalent  to  the  first  order  basis  functions  of  Wang  and  Webb  [14], 
but  are  presented  here  in  a  simpler  formulation  that  allows  existing 
RWG  codes  to  be  easily  modified  to  use  them. 

2.  DESCRIPTION  OF  THE  PROBLEM 

Whra  we  discretize  a  dosed  surface  using  Nf  triangular  patrheB,  we 
obtain  a  set  of  =  2Nf/2  edges  connecting  these  triangles.  We  want 
to  develop  a  set  of  basis  functions  that  can  represent  an  arbitrary  linear 
current  distribution  over  each  triangle.  For  eadi  triangle  »  we  would 
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Figure  1.  Linear  current  components  of  two  adjacent  triangles. 


have  a  different  current  distribution  given  by: 

Vi)  —  {p-xi^i  +  ^xiVi  +  Pei)  +  {Pyi^i  +  +  Cyj)  yi  (1) 

where  Xi  and  j/j  are  local  coordinates  in  the  directions  tangential  to 
the  triangle  i.  Since  we  have  Nf  trieingles,  the  number  of  degrees  of 
freedom  for  this  problem  is  equal  to  6fVy. 

For  electromagnetic  problems  we  are  interested  only  in  current 
distributions  whose  components  normd  to  the  edges  are  continuous, 
to  avoid  the  presence  of  artificial  line  charges.  From  Fig.  1,  we  see 
that  by  imposing  this  additional  condition  for  each  edge  of  our  surface 
we  are  adding  two  constraints  for  each  edge  {oni  =  Onj,  Cni  =  nnj)t 
thereby  reducing  the  degrees  of  fireedom  firom  6Nf  to  SNf.  We  ne^, 
consequently,  a  set  of  3JVy  divergence-conforming  basis  functions  [15] 
to  properly  solve  this  problem. 

In  [1],  Rao  et  al  propose  a  set  of  basis  functions  (RWG  basis) 
having  one  function  associated  with  each  edge,  leading  therefore  to  a 
set  of  only  ZNf/2  basis  functions.  Clearly,  su^  set  cannot  represent 
ein  arbitrary  linear  current  distribution  as  described  earlier,  i.e.,  it  is 
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Figure  2.  Triangle  pair  and  geometrical  parameters  associated. 


not  complete  in  that  sense.  This  has  already  been  pointed  out  in  [14] 

and  [15].  r  i-  u  • 

To  solve  this  problem  we  developed  a  new  set  of  Imear  basis 

functions  (which  we  shall  call  the  fet  order  basis)  where  we  have  two 
vector  basis  functions  associated  with  each  edge  b-c,  as  seen  in  Fig.  2. 


/i(0  = 


I 


(2A+)2 

(2A-)2 

I 

(2A+)2 


1^ 


f  inT^ 
f  inT" 

r  in 
r  in  T~ 


(2) 


where  p"*"  =  f  —  ro+,  p  —  —  f  and  p^J.  —  fa*---  ^ 

the  area  of  the  respective  triangle  and  “x”  denotes  the  cross  product. 
Each  of  these  basis  functions  is  parallel  to  one  external  edge  and  has  a 
magnitude  linearly  proportional  to  the  distance  to  the  other  external 
edge,  as  shown  in  Fig.  3. 

The  first  order  basis  functions  share  the  following  properties  with 
the  RWG  functions: 

•  The  current  has  no  component  normal  to  the  boundary  of  the 
surface  formed  by  triangles  T+  and  T~  (external  edges). 

•  The  component  of  current  normal  to  the  edge  between  the  two 
triangles  (internal  edge)  is  continuous  across  the  edge.  While  it 
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current 

intensity 


Figure  3.  First  order  trianguleir  patch  basis  function  (of  type  1) 
intensity  and  direction. 


was  a  constant  for  the  RWG  basis  functions,  it  is  now  a  linear 
function,  going  from  zero  at  one  vertex  to  one  at  the  opposite 
vertex. 

•  The  surface  divergence  of  these  functions,  which  is  proportional 
to  the  surface  charge  density,  is  constant  over  each  triangle.  The 
resulting  surface  charge  density  has  the  same  magnitude  but 
opposite  signs  on  the  two  triangles. 

In  fact,  it  can  be  shown  that  an  RWG  basis  function  can  be 
obtained  by  adding  the  two  first  order  basis  functions  of  equal  strength 
for  the  same  edge.  Because  of  the  properties  described,  it  will  be  shown 
that  the  electric  field  integral  equation  (EFIE)  formulation  used  wth 
the  RWG  functions  [1]  can  still  be  used,  with  only  minor  modifications 

for  the  first  order  basis  functions. 

As  an  example  of  the  improvement  that  can  be  obtained  with  the 
first  order  basis  functions.  Fig.  4  shows  the  best  representation  (in  the 
mean  square  sense)  of  a  sinusoidal  cmrrent  on  a  rectangular  plate  that 
can  be  obtained  with  the  RWG  and  the  first  order  basis  functions  with 
the  gnTnp.  number  of  triangular  patches.  We  can  see  that  the  first  order 
basis  representation  is  significantly  smoother  and  closer  to  the  exact 
expression  than  the  RWG  basis. 


Trintmalia-Ling 

Rao-Wilton-Glisson 


First  order  triangular  patch  basis  functions 


1527 


0.5  -05 


(c) 

Figure  4.  Sinusoidal  current  representation  using  first  order  basis 
functions  (b)  and  the  original  RWG  basis  functions  (c)  over  a  square 
plate,  using  72  triangular  patches.  Figure  (a)  shows  the  A- A'  cut  of 
both  representations. 


3.  EFEE  FORMULATION 

We  will  now  outline  the  implementation  of  the  EFIE  using  the  first 
order  basis.  Using  the  standard  EFIE  formulation  with  the  Galerkin 
method,  we  obtain: 

{*^5/7711,2^  ~  /i7ii,2^  “I"  /mi, 2^  (3) 

where  (a, M  ^  f  a^bdS  and  fmi  2  denotes  the  basis  functions  in  (2)  at 

the  m-th  edge.  As  in  [1],  because  of  the  properties  of  fm  at  the  edges, 
the  last  term  in  (3)  can  be  rewritten  as 

(V$,  frma)  =  V  •  fX,,  I  ^dS  +  V-  j  ^dS  (4) 


'  fm  ~  i 
•/mi  2 


2A±' 


with 


(5) 


1528 


TVintinalia  and  Ling 


And 

=  j  A  -  fmi,2d'S  +  j  A  -  fmi,2dS,  (6a) 

r+ 

fmi,^  =  J  ^ j  ^ 

r+  r;; 

One  difference  should  be  noted  here.  While  in  [1]  the  values  of  the 
t6rins  involving  the  vector  potential  and  the  incident  field  in  (6)  were 
approximated  by  their  values  at  the  centroid  of  the  triangles,  with  our 
new  set  of  basis/testing  functions  one  should  use  a  different  and  more 


precise  approximation: 

fA  frm,2dS  ^ 

(7a) 

fp± 

■*m 

^  i=l 

f^-U,2dS  S 

/Tl± 

i=l 

(7b) 

with 

ra  =  ^ra  +  -n  +  -rc, 

1-  1-  2^ 
ra  =  gr,  +  -r6  +  -rc. 


(8) 


This  is  necessary  because  we  have  3  degrees  of  freedom  for  each  triangle 
now,  so  we  need  at  least  3  testing  points.  Moreover,  this  approximation 
leads  to  exact  results  when  the  vector  potential  or  the  incident  field 
varies  linearly.  Similarly,  the  scalar  potential  integral  in  (4)  can  be 


approximated  as: 


u 


(9) 


The  remaining  steps  to  derive  the  matrix  equation  are  similar 
to  those  in  [1].  In  particular,  the  integrals  related  to  the  scalar  and 
vector  potential  determination,  after  transforming  to  a  local  system 
of  normalized  area  coordinates,  can  be  written  as  a  functions  of  the 
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integrak  if,  i?«  and  if: 

"1.2  jAirue  ’ 


(10) 


= 

■^i.a 


where 

pQ 


rP9 

R 

Here,  the  notation  ^  or  >1^  ^  means  the  potential  at  a  point  r  in 
the  triangle  p  due  to  basis  function  (1  or  2)  associated  with  edge  n 
flowing  on  triangle  q.  The  positions  fi,,  f2q,  fzq  are  the  three  vertices 
of  triangle  q  (independent  of  n)  eind  the  other  parameters  are  those 
specified  in  Fig.  2. 

4.  NUMERICAL  RESULTS 

In  order  to  compare  the  convergence  of  the  first  order  basis  against  the 
standard  RWG  basis  we  first  analyzed  the  scattering  by  a  perfectly 
conducting  square  plate.  Fig.  5  and  Fig.  6  show  respectively  the 
predicted  co-polarization  and  cross-polarization  RCS  of  a  square  plate 
of  side  lA  at  normal  incidence,  plotted  as  a  function  of  the  number  of 
basis  functions  used  {ZNf/2  for  RWG  and  ZNf  for  first  order  basis). 
We  can  see  a  much  faster  convergence  with  the  first  order  basis  than 
with  the  RWG  basis.  Note  that  this  is  so  even  though  for  the  same 
number  of  bases,  the  triangul2U'  mesh  for  the  RWG  result  is  twice  as 
dense  as  that  used  in  the  first  order  basis  result. 

Fig.  7  shows  the  current  distribution  on  the  plate  obtained  using  a 
discretization  of  72  triangular  patches  for  both  sets  of  basis  functions. 
We  observe  rather  abrupt  variations  of  the  current  along  the  transversal 
direction  for  the  RWG  basis.  For  the  first  order  basis,  on  the  other 


[ll(f„  -  f.±)  X  ft, ,11  Jf  +  ||(r„  -  X 
+  (11) 


1  l-»7  ,t.p 

I  I  R 

1  l-»7  -up 

r  t  p-jKH 

h  R 

0  0 

0  0 

]  h  R 

n  n 

=  (12) 

u  u 

|r-r'|,  r'€r«. 

feT^. 
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Fieure  5.  Predicted  RCS  for  normal  incidence  at  a  square  plate  of 
side  1  m  (A  =  Im),  as  a  function  of  the  number  of  basis  functions 

used. 


Fieure  6.  Predicted  cross-polarization  for  normal  incidence  at  a 
square  plate  of  side  lA,  as  a  function  of  the  number  of  basis  functions 

used. 
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(b) 

Figure  7.  Magnitude  of  current  distribution  (x  component)  over  a 
PEC  square  plate  of  side  lA  obtained  using  a  72  triangular  pat<^ 
discretization  for  (a)  RWG  basis  functions  and  (b)  first  order  basis 
functions. 

hand,  the  current  behavior  is  smoother,  as  expected.  Also,  for  the  &st 
order  basis  solution  the  current  distribution  reaches  higher  values  close 
to  the  edges,  as  it  should  due  to  the  edge  singularity. 

As  a  second  example  we  examined  the  scattering  from  a  sphere. 
Fig.  8  shows  the  current  distribution  over  a  PEC  sphere  u^g  a 
discretization  of  160  triangular  patches.  In  this  figure  only  tbe  r^ 
part  of  the  B  component  is  shown,  but  similar  behavior  is  exhibited 
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Figure  8.  Current  distribution  (real  part  only)  over  a  PEC  sphere 
of  radius  0.4  A,  obt^dned  using  a  160  triangular  patch  discretization 
for  RWG  bAsis  and  first  order  basis  (TL),  compared  against  the  exact 
Mie’s  solution.  The  wave  is  incident  firom  9 —  OP. 


0  03  1  13  ?  23  3  3.5  *  AS  5 

*10"^ 


Figure  9.  Real  part  of  the  total  current  density  (vector  magnitude) 
over  a  PE5C  sphere  of  radius  0.4  A,  obtained  using  a  160  triangular 
patch  discretization  for  RWG  basis  functions  (a)  and  first  order  basis 
functions  (b).  The  wave  is  incident  firom  0  =  0“  (top). 
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Fierure  10.  Mean  square  error  in  current  distribution  for  a  0.4  A  radius 
sphere,  obtained  using  the  RWG  basis  and  the  first  order  basis,  as  a 
function  of  the  number  of  triangular  patches  used. 

for  the  imaginary  part  as  weU  as  the  other  current  component.  We 
can  see  that  the  agreement  with  the  exact  Mie’s  solution  is  better  for 
the  first  order  basis  than  for  the  RWG  basis.  Fig.  9  shows  the  r^ 
part  of  the  total  current  density  on  the  surface  of  the  sphere.  Notice 
that  the  first  order  basis  solution  presents  a  much  smoother  behavior 
than  the  RWG  one.  Fig.  10  shows  the  total  mean  square  error  m  the 
current  distribution  obtained  using  both  formulations,  as  a  function  of 
the  number  of  facets,  for  the  same  sphere  problem.  Here  we  observe 
^gAin  tbftt  the  convergence  of  the  first  order  basis  solution  is  much 

faster  than  the  RWG  solution.  ^  ^ 

We  also  ran  a  frequency  scan  of  the  RCS  prediction  for  a  0.2  m 
radius  sphere  using  both  bases.  The  results,  compared  against  the 
exact  Mie’s  solution,  are  shown  in  Fig.  11.  We  note  that  at  lower 
frequencies  (discretization  <  0.2  A)  both  models  provide  very  good 
agreement  with  the  exact  solution.  For  higher  frequencies  though 
Id  >  0.3  A)  both  models  start  to  diverge  from  the  ocact  solution,  but 
the  first  order  basis  is  more  stable  than  the  RWG  basis  as  the  frequency 
increases. 
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Figure  11.  RCS  prediction  as  a  function  of  frequency  for  a  sphere 
of  radius  0.2  m,  using  a  424-facet  model  (discretization  S  0.05  m)  for 
RWG  basis  fimctions  and  first  order  basis  (TL)  functions,  compared 
against  the  exact  Mie’s  solution. 

5.  SUMMARY 

In  this  paper  we  presented  an  enhanced  set  of  triangular  patch  basis 
functions  which  represent  an  improvement  over  the  traditional  Rao, 
Wilton  and  Glisson  triangular  patch  basis  functions.  These  basis 
functions  axe  equivalent  to  the  first  order  set  described  in  [14]  but  are 
presented  here  in  a  simpler  formulation  that  allows  its  use  in  existing 
RWG  codes  with  minimal  modifications.  It  has  been  shown  that 
these  new  basis  functions  provide  a  much  better  representation  of  the 
current  distribution  than  RWG  while  keeping  the  same  computational 
complexity.  In  the  examples  we  have  analyzed,  we  always  obtained, 
for  the  same  triangiilation  scheme,  a  better  result  using  the  first 
order  basis  functions  in  terms  of  RCS  prediction.  In  terms  of  current 
distribution,  the  comparative  performance  is  even  better  for  this  new 
formulation.  Therefore  we  can  conclude  that  this  set  of  basis  functions 
really  represents  an  improvement  over  the  traditional  RWG  basis. 
Also,  its  implementation  can  be  easily  extended  to  the  magnetic  field 
and  combined  field  integral  equation  formulations. 

As  a  final  remark  we  should  mention  that  the  development  of 
hi^er-order  basis  functions  has  received  much  attention  recently  [14- 
l"^.  Higher-order  basis  achieves  a  faster  convergence  rate  at  the 
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price  of  increased  computational  cost  for  the  matrix  elements.  In 
comparison,  the  formulation  presents  here  does  not  claim  to  provide 
faster  convergence  rate,  but  is  a  fast  and  easy  modification  to  existing 
RWG-based  computer  codes  to  achieve  enhanced  performance. 
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Abstract 

The  active  element  patterns  of  antenna  arrays  are  strongly  affected  by  the  interaction 
between  the  elements  and  the  mounting  platform.  In  this  paper,  a  combined  model  is 
proposed  to  model  the  array  response  including  both  mutual  coupling  and  the  platform 
effect.  The  mutual  coupling  is  described  by  a  coupling  matrix  while  the  platform 
scattering  is  represented  using  the  point  scatterer  model.  An  algorithm  based  on  the 
matching  pursuit  concept  is  proposed  to  determine  the  model  coefficients.  Numerical 
results  show  that  the  active  element  patterns  are  well  described  by  the  model. 


Key  words:  Array  mutual  coupling,  array-platform  interaction,  point  scatterer  model 


L  Introduction 


It  is  well  known  that  the  chciracteristics  of  an  antenna  array  are  strongly  affected  by 
the  mutual  coupling  between  array  elements  [1-8].  Mutual  coupling  is  often  modeled 
using  a  coupling  matrix  [3-5].  For  an  AT-element  array,  the  coupling  matrix  is  an  by 
matrix  relating  the  active  element  patterns,  which  are  the  antenna  element  patterns  in  the 
array  environment,  to  the  free-standing  element  patterns.  Under  most  conditions,  the 
coupling  matrix  provides  a  very  effective  way  to  model  the  array  response  in  the  presence 
of  mutual  coupling.  However,  if  the  array  is  mounted  on  a  platform,  the  active  element 
patterns  are  further  modified  by  the  scattering  from  the  platform  [9].  In  this  case,  the 
mutual  coupling  matrix  alone  is  insufficient  to  describe  the  radiation  physics.  Recently, 
we  have  shown  that  antenna-platform  interactions  can  be  well  described  by  a  sparse  point 
scatterer  model  [10,11].  In  this  approach,  the  scattered  field  from  the  platform  is 
approximated  by  the  contributions  of  a  few  point  scatterers.  In  this  paper,  we  propose  a 
combined  model  that  takes  into  account  of  both  the  mutual  coupling  and  the  antenna- 
platform  interaction.  This  model  combines  the  coupling  matrix  approach  and  the  point 
scatterer  description  of  the  platform  scattering  to  fully  characterize  the  radiation  physics 
associated  with  arrays  mounted  on  complex  platforms.  In  Section  H,  the  mutual  coupling 
model  is  reviewed  and  the  combined  model  is  described.  In  Section  HI,  an  algorithm 
based  on  the  matching  pursuit  concept  [12]  is  proposed  to  determine  the  model 
parameters  including  the  position  and  strength  of  the  point  scatterers.  Numerical  results 
are  provided  in  Section  IV  to  demonstrate  the  effectiveness  of  the  model. 


11.  The  Combined  Model 


Assume  an  array  with  N  elements  is  located  in  the  Ay-plane  and  receives  a  plane  wave 
incident  from  angle  ^  in  the  azimuth  plane.  The  standard  mutual  coupling  model  [3-5] 
relates  the  active  element  patterns  to  the  free-standing  element  patterns  by  a  coupling 
matrix  as  follows: 


A 


true 


=  CA 


theo 


(1) 


where  ktrue  and  \.,heo  are  the  active  and  free-standing  element  patterns,  respectively,  and 
C  is  the  Nhy  N  coupling  matrix.  For  antenna  elements  having  isotropic  patterns  in  the  xy- 
plane,  the  elements  of  ktheo  are  given  by: 


(2) 


where  (x,,  y,)  is  the  coordinate  of  the  ith  element.  The  active  pattern  matrix  has  the 
same  structure  as  k,heo,  except  each  row  represents  the  antenna  element  pattern  in  the 
array  environment.  Given  the  active  element  patterns  at  a  few  incident  angles,  the  matrix 
C  can  be  determined  using  the  pseudo  inverse  as  follows: 

c  =  a„a“,(a,^.a1)-'  (3) 


In  [8],  we  demonstrated  that  the  model  in  (1)  is  valid  if  the  current  distribution  on 
each  stand-alone  element  has  a  constant  shape  for  all  incident  angles.  It  was  also  shown 
that  when  the  antenna  elements  are  resonant,  the  constant-current-shape  condition  is 
approximately  satisfied.  This  situation  should  still  hold  true  even  if  a  platform  exists. 
Thus  the  mutual  coupling  model  can  be  extended  to  the  case  when  a  platform  is  present, 
provided  that  we  replace  the  free-standing  element  patterns  in  ktheo  by  the  array  element 
patterns  in  the  presence  of  the  platform  (but  in  the  absence  of  the  other  elements).  Next 
we  consider  how  the  platform  affects  the  stand-alone  pattern  of  each  element. 


For  a  stand-alone  element,  the  platform  scattering  changes  the  incident  field  over  the 
element.  The  total  incident  field  on  the  element  can  be  viewed  as  a  sum  of  the  direct 
incident  field  and  the  scattered  field  from  the  platform.  As  shown  in  [10,11],  for  an 
electrically  large  platform,  the  scattered  field  can  be  adequately  described  using  a  sparse 
point  scatterer  model.  Consequently,  we  can  incorporate  the  platform  effect  in  the  mutual 
coupling  model  as  follows: 


^true  ~  ^^^theo  "*■  plat  ) 


(4) 


In  this  model,  Apiat  represents  the  incident  fields  upon  the  point  scatterers  representing  the 
platform.  It  has  the  exact  same  form  as  an  ideal  element  pattern  in  (2),  where  (x,,  y,) 
represents  the  coordinate  of  the  point  scatterer.  Each  column  of  the  matrix  P  describes  the 
scattering  from  each  point  scatterer  to  the  array  elements.  It  includes  both  the  strength  of 
the  scatterer  and  the  phase  delay  from  the  scatterer  to  the  array  elements.  I  is  the  identity 
matrix.  It  is  clearly  shown  in  this  equation  that  the  stand-alone  element  patterns  consist  of 
two  parts,  the  free-standing  element  patterns  and  the  modification  due  to  the  platform. 
Thus  the  model  is  consistent  with  the  physics  of  the  actual  problem. 


III.  Determination  of  Model  Coefficients 

To  determine  the  positions  of  the  point  scatterers,  one  standard  way  is  to  generate  a 
projection  image  [13]  in  two-dimensional  (2D)  space.  In  this  method,  we  assume  a  point 
scatterer  is  located  at  (x,  y).  The  matrix  Apiat  becomes  a  row  vector  given  by: 


jk(xcos^+y&m^) 


(5) 


We  can  then  determine  the  matrix  C[I  P]  using  the  pseudo  inverse  technique  similar  to 
(3).  Once  C  and  CP  are  explicitly  obtained,  the  column  vector  P  can  be  determined  by 
left-multiplying  CP  by  C‘^.  The  above  process  can  be  viewed  as  a  projection  of  the 
assumed  point  scatterer  response  onto  the  actual  array  response  including  platform 
effects.  The  energy  in  P  will  be  maximum  if  we  choose  the  correct  point  scatterer 
location.  Thus  we  can  take  the  norm  of  P  to  form  a  2D  image  as  a  function  of  (x:,  y).  The 
positions  of  the  point  scatterers  are  determined  from  the  peaks  in  the  image. 

However,  the  drawback  of  this  approach  is  that  the  2D  image  is  generated  using  one¬ 
dimensional  data.  The  sidelobes  of  the  peaks  are  expected  to  be  very  high  such  that  the 
peaks  corresponding  to  the  weaker  point  scatterers  may  be  shadowed.  To  illustrate  this, 
we  consider  the  simplified  projection  image  of  a  single  point  scatterer  without  any  array 
elements.  Suppose  the  point  scatterer  is  located  at  (xo,  yo),  the  image  strength  at  an 
assumed  point  (jc,  y)  is  given  by 


jk{xcos<p+ysin  ^(J^ocos^+yosin  (ft) 

jkr{cQ$  ipQ  cos  ^+sin  <pQ  sin 


=  \27a,{kr)\ 


(6) 


where  r  =  the  distance  between  (xr,  y)  and  (xo,  yo)<  Thus  the 

projection  image  of  a  single  point  scatterer  is  a  zeroth  order  Bessel  function  centered  at 
ixo,yo).  The  sidelobe  of  the  Bessel  function  is  so  high  (about  8  dB)  that  simultaneous 
identification  of  point  scatterers  is  nearly  impossible.  To  overcome  this  problem,  we 
employ  a  matching  pursuit  algorithm  [12]  to  determine  the  point  scatterer  positions 
iteratively. 


In  the  matching  pursuit  algorithm,  a  projection  image  is  generated  at  each  iteration. 
The  first  image  is  generated  using  IIPII  as  described  above.  The  strongest  point  in  the 
image  is  picked  to  be  the  first  point  scatterer.  The  position  information  is  stored  in  Aput. 
In  the  second  iteration,  we  assume  the  position  of  the  next  point  scatterer  and  append  Apiat 
by  one  row.  A  new  projection  image  is  then  formed  using  the  norm  of  the  last  column  of 
P  and  the  second  point  scatterer  is  determined.  The  process  is  repeated  until  the  total 
energy  in  the  projection  image  is  small  enough  to  be  neglected. 

Using  this  method,  the  point  scatterers  are  determined  one  at  a  time.  Once  a  point 
scatterer  is  included  in  the  model,  its  sidelobe  is  taken  into  account  so  that  it  will  not 
overshadow  the  weaker  point  scatterers  in  the  subsequent  projection  images. 

IV.  Numerical  Results 

In  this  section,  two  examples  are  provided  to  validate  the  combined  model  of  mutual 
coupling  and  platform  effects.  In  the  first  example,  we  simulate  a  linear  array  of  half¬ 
wave  dipoles  mounted  on  a  ship-like  platform.  The  array  and  the  platform  geometry  are 
shown  in  Fig.  1.  The  separation  between  the  array  elements  is  0.45  wavelength.  We 
simulate  the  problem  with  FISC  [14],  which  is  a  method  of  moments  code  based  on  the 
multi-level  fast  multipole  method.  The  active  element  patterns  at  36  uniformly  distributed 
incident  angles  are  used  to  extract  the  model  coefficients  in  (4).  Fig.  2  shows  the 
projection  image  generated  in  the  first  iteration  of  the  matching  pursuit  algorithm.  A 
strong  point  scatterer  can  be  identified  at  the  back  comer  of  the  front  cylinder.  The  high 
sidelobes  associated  with  this  scatterer  can  also  be  seen  from  the  image.  In  Fig.  3,  the 
first  15  point  scatterers  extracted  by  the  matching  pursuit  algorithm  are  plotted  over  the 


platform.  The  scatterers  are  mostly  distributed  along  the  centerline  of  the  ship,  and 
correspond  to  the  expected  specular  reflection,  tip  diffraction  and  comer  scattering 
mechanisms  on  the  topside  stmcture.  Two  weak  scatterers  are  asymmetric  about  the 
center  line.  The  asymmetry  is  due  to  numerical  noise  in  the  extraction. 

Once  the  combined  model  is  obtained,  we  can  reconstmct  the  active  element  patterns 
at  arbitrary  incident  angles.  The  reconstructed  and  simulated  active  element  patterns  are 
plotted  in  Fig.  4,  in  dashed  and  solid  lines,  respectively.  The  agreement  between  the  two 
is  very  good.  As  a  comparison,  we  also  obtained  a  model  based  on  the  traditional  mutual 
coupling  model  in  (1)  without  any  platform  considerations.  The  active  element  patterns 
reconstmcted  from  the  traditional  model  are  plotted  in  dash-dot  lines.  It  can  be  seen  that 
the  fine  stmctures  of  the  array  response  are  not  captured  by  the  traditional  model.  The 
difference  is  more  dramatic  when  we  examine  the  phase  error  of  the  reconstructed 
element  patterns  by  the  combined  model  and  the  traditional  model.  They  are  shown  in 
Fig.  5  in  solid  and  dashed  curves,  respectively.  We  can  see  that  the  combined  model 
performs  considerably  better  in  terms  of  phase  error. 

In  some  applications  such  as  direction  finding,  the  phase  of  the  active  element 
patterns  can  be  quite  important.  To  demonstrate  this,  we  simulate  the  array  response 
when  there  are  two  incoming  signals  separated  by  15°.  A  signal-to-noise  ratio  of  30  dB  is 
assumed.  The  MUSIC  superresolution  algorithm  is  employed  to  determine  the  incoming 
signal  directions  in  the  presence  of  noise.  The  normalized  MUSIC  power  spectrum  using 
the  exact  active  element  patterns  is  plotted  as  solid  curve  in  Fig.  6.  The  power  spectra 
resulting  from  the  combined  model  and  the  traditional  coupling  model  are  plotted  in 


dashed  and  dash-dot  curves,  respectively.  We  can  see  that  the  two  incident  signals  can  be 
well  separated  using  the  combined  model,  while  the  traditional  model  fails. 

In  the  second  example,  a  7-element  circular  array  of  dipoles  with  a  center  pole  and 
some  surrounding  structures  is  simulated.  The  array-platform  geometry  is  shown  in  Fig. 
7.  The  combined  model  is  extracted  using  the  active  element  patterns  computed  at  36 
angles  and  is  used  to  reconstruct  the  patterns.  The  projection  image  in  the  first  iteration  of 
the  matching  pursuit  algorithm  is  shown  in  Fig.  8.  The  strongest  scattering  from  the 
platform  is  observed  to  be  coming  from  the  comer  of  the  nearby  box.  The  extracted  point 
scatterers  are  also  plotted  as  circles  in  this  figure.  While  the  ring  stmcture  in  the  image 
resembles  the  case  where  there  is  a  point  scatterer  at  the  center  of  the  array,  none  is 
extracted  by  the  algorithm.  We  believe  this  is  because  the  scatterers  are  too  close  to  the 
antenna  elements.  As  a  result,  Atheo  and  Apiat  in  (4)  are  too  similar  and  the  effect  from  the 
point  scatterers  is  included  in  the  C  matrix  instead.  Despite  this  degeneracy,  the  resulting 
model  is  still  adequate  in  describing  the  active  element  patterns.  Fig.  9  shows  the  phase 
errors  in  the  reconstmcted  active  element  patterns  based  on  the  model.  We  again  see  that 
the  combined  model  performs  much  better  than  the  traditional  coupling  model.  The 
direction  finding  results  corresponding  to  the  different  active  element  patterns  are  plotted 
in  Fig.  10.  As  expected,  the  combined  model  outperforms  the  traditional  coupling  model, 
as  it  can  clearly  distinguish  the  two  incoming  signals. 

V.  Conclusions 

In  this  paper,  we  have  proposed  a  combined  model  to  describe  the  response  of  an 
antenna  array  in  the  presence  of  both  mutual  coupling  and  platform  effects.  The  mutual 


coupling  between  the  array  elements  is  described  by  a  coupling  matrix.  In  addition,  the 
platform  effect  is  approximated  using  the  point  scatterer  model.  We  have  also  proposed 
an  algorithm  based  on  matching  pursuit  to  determine  the  model  coefficients.  Numerical 
results  show  that  the  active  element  patterns  reconstructed  from  the  model  is  in  good 
agreement  with  the  actual  data.  The  main  attractiveness  of  the  new  model  is  that  it 
provides  physical  insights  into  the  coupling  and  scattering  mechanisms.  By 
parameterizing  the  active  element  patterns  using  the  proposed  model,  we  can  locate  the 
position  and  strength  of  the  dominant  point  scatterers  that  give  rise  to  platform  scattering. 

Acknowledgment 

This  work  was  supported  by  the  Office  of  Naval  Research  under  Contract  No. 


N00014-98-1-0178  and  NOOO 14-0 1-0234. 


References 


[1]  D.  F.  Kelley  and  W.  L.  Stutzman,  “Array  antenna  pattern  modeling  methods  that 
include  mutual  coupling  effects,”  IEEE  Trans.  Antennas  Propagat.,  v.  41,  pp  1625- 
1632,  Dec.  1993. 

[2]  I.  J.  Gupta  and  A.  A.  Ksienski,  “Effect  of  mutual  coupling  on  the  performance  of 
adaptive  arrays,”  IEEE  Trans.  Antennas  Propagat.,  v.  31,  pp  785-791,  Sept.  1983. 

[3]  B.  Friedlander  and  J.  Weiss,  “Direction  finding  in  the  presence  of  mutual  coupling,” 
IEEE  Trans.  Antennas  Propagat.,  v.  39,  pp  273-284,  Mar.  1991. 

[4]  T.  Svantesson,  “Modeling  and  estimation  of  mutual  coupling  in  a  uniform  linear 
array  of  dipoles,”  IEEE  Int.  Conf.  Acoust.  Speech  Signal  Processing,  pp.  2961- 
2964,  Phoenix,  1999. 

[5]  C.  M.  S.  See,  “A  method  for  array  calibration  in  parametric  sensor  array 
processing,”  IEEE  Singapore  Int.  Conf.  on  Communication  Systems,  pp.  915-919, 
Singapore,  1994. 

[6]  R.  S.  Adve  and  T.  K.  Sarkar,  “Compensation  for  the  effects  of  mutual  coupling  on 
direct  data  domain  adaptive  algorithms,”  IEEE  Trans.  Antennas  Propagat.,  vol.  48, 
pp.  86-94,  Jan.  2000. 

[7]  T.  Su,  K.  Dandekar  and  H.  Ling,  “Simulation  of  mutual  coupling  effect  in  circular 
arrays  for  direction  finding  applications,”  Microwave  Optical  Tech.  Lett.,  vol.  26, 
pp.  331-336,  Sept.  2000. 

[8]  T.  Su  and  H.  Ling,  “On  modeling  mutual  coupling  in  antenna  arrays  using  the 
coupling  matrix,”  Microwave  Optical  Tech.  Lett.,  vol.  28,  pp.  231-237,  Feb.  2001. 


[9]  T.  Milligan,  F.  Obelleiro,  L.  Landesa,  J.  Taboada  and  J.  Rodriguez,  “S)nithesis  of 
onboard  array  antennas  including  interaction  with  the  mounting  platform  and 
mutual  coupling  effects,  ”  IEEE  Antennas  Propagat.  Mag.,  vol.  43,  pp.  76-82,  Apr. 
2001. 

[10]  C.  Ozdemir,  R.  Bhalla  and  H.  Ling,  “A  radiation  center  representation  of  antenna 
radiation  pattern  on  a  complex  platform,”  IEEE  Trans.  Antennas  Propagat.,  vol. 
AP-48,  pp.  992-1000,  June  2000. 

[11]  T.  Su,  C.  Ozdemir  and  H.  Ling,  “On  extracting  the  radiation  center  representation 
of  antenna  radiation  patterns  on  complex  platform,”  Microwave  Optical  Tech.  Lett, 
vol.  26,  pp.  4-7,  July  2000. 

[12]  S.  G.  Mallat  and  Z.  Zhang,  “Matching  pursuits  with  time-frequency  dictionaries,” 
IEEE  Trans.  Signal  Processing,  v.  41,  pp.  3397-3415,  Dec.  1993. 

[13]  G.  T.  Herman,  Image  Reconstruction  from  Projections,  Academic  Press,  New  York, 
1980. 

[14]  User’s  Manual  for  FISC  {Fast  Illinois  Solver  Code),  Center  for  Computational 
Electromagnetics,  University  of  Illinois  at  Urbana-Champaign  and  DEMACO,  Inc., 


IL,Jan.  1997. 


Fig.  1.  Array-platform  geometry  of  example 
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Fig.  4.  Actual  and  reconstructed  active  element  patterns 
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- Combined  model 

- Traditional  coupling  model 
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Fig.  6.  MUSIC  power  spectra  using  differ 

- Actual  data 

- Combined  model 

- Traditional  coupling  mode 


Fig.  8.  First-iteration  projection  image  and  extracted  point  scatterers 
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(b)  element  2 


Fig.  9.  Phase  error  of  the  reconstructed  active  element  patterns 
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Fig.  10.  MUSIC  power  spectra  using  different  active  element  patterns 
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Abstract: 

This  paper  investigates  the  benefit  of  mutual  coupling  compensation  via  a  Method 
of  Moments  (MoM)  approach  in  a  uniform  circular  antenna  array  operating  at  1 .8  GHz. 
This  mutual  coupling  compensation  technique  is  applied  to  a  direction  of  arrival  (DOA) 
study  of  up  to  two  co-channel  mobile  users.  Field  measurements  and  computer 
simulations  are  examined  to  explore  the  assumptions  of  the  technique  and  verify  its  effect 
when  using  the  Bartlett  and  MUSIC  DOA  algorithms.  Computer  simulations  considering 
the  application  of  the  technique  to  downlink  beamforming  are  also  included. 
Experimental  results  show  that  the  mutual  coupling  compensation  technique  improves 
uplink  DOA  algorithm  performance  primarily  by  reducing  unwanted  sidelobe  levels. 
This  reduction  in  sidelobe  levels  aids  in  downlink  beamforming  weight  design. 
Specifically,  simulation  results  show  that  use  of  the  compensation  technique  allows 
DOA-based  downlink  beamforming  algorithms  to  perform  similarly  to  spatial  signature- 
based  algorithms.  All  field  measurements  were  made  using  the  smart  antenna  testbed  at 
the  University  of  Texas  at  Austin. 

Index  Terms:  Adaptive  Arrays,  Antenna  Array  Mutual  Coupling,  Array  Signal 
Processing,  Direction  of  Arrival  Estimation 


1.  Introduction: 


Planning  for  3"^^  generation  cellular  systems  is  currently  underway  to  handle  the 
increasing  demand  for  wireless  services.  Smart  antenna  technology,  making  use  of 
adaptive  antenna  arrays,  is  capable  of  improving  system  capacity  and  quality,  and  will  be 
an  essential  component  in  these  next  generation  cellular  systems.  Fully  realizing  the 
potential  of  smart  antenna  systems  requires  an  appreciation  of  the  underlying 
electromagnetic  effects  governing  antenna  array  operations.  These  electromagnetic 
effects  include  mutual  coupling  as  well  as  amplitude  and  phase  mismatch  between  array 
elements. 

It  has  become  clear  to  the  smart  antenna  community  that  the  actual  response  of 
the  antenna  array  can  deviate  significantly  from  the  assumed  and  simplistic  model  due  to 
electromagnetic  coupling  among  the  antenna  elements  as  well  as  scattering  from  the 
antenna  tower  and  nearby  structures.  Most  studies  of  mutual  coupling  in  antenna  arrays 
can  be  best  characterized  as  array  calibration  methods  [1-5].  These  methods  are 
generally  derived  using  parametric  antenna  array  models  that  do  not  directly  consider  the 
underlying  electromagnetic  principles.  A  costly  alternative  to  these  methods  involves 
determining  the  actual  array  response  using  field  measurements. 

The  same  problem  can  also  be  approached  from  a  fundamental  electromagnetics 
perspective  [6-10].  In  this  paper,  this  approach  is  applied  to  uniform  circular  arrays  to 
demonstrate  in  field  measurements  the  effects  of  mutual  coupling  compensation.  In  this 
mutual  coupling  compensation  technique,  the  Method  of  Moments  [11]  (MoM)  is  used  to 
compute  an  array  response  that  compensates  for  mutual  coupling.  The  paper  will  then 
examine  uplink  (mobile  to  basestation)  DOA  algorithm  performance  in  field 


measurements  for  scenarios  involving  up  to  two  co-channel  mobile  users.  The  study  will 
also  consider  compensation  in  downlink  (basestation  to  mobile)  beamforming  scenarios 
and  thus  determine  the  applicability  of  the  technique. 

Section  n  introduces  a  mathematical  model  of  the  smart  antenna  array  and 
describes  the  manner  in  which  the  mutual  coupling  compensation  technique  can  be 
derived.  Section  in  briefly  describes  the  MoM  calculation  necessary  for  coupling 
characterization.  Field  measurement  and  simulation  results  of  the  considered 
compensation  technique  are  discussed  in  Section  IV. 

11.  Smart  Antenna  Model 
A.  Array  Output 

Assume  that  there  are  M  elements  in  the  base  station  antenna  array  located  at  (Xi, 

yi),  1  <  i  <  M.  We  define  the  concept  of  a  steering  vector  (sometimes  called  array 

response  vector)  that  can  be  thought  of  as  the  spatial  analog  of  an  impulse  response  in 

temporal  processing.  Specifically,  it  characterizes  the  relative  phase  response  of  each 

antenna  array  element  to  an  incident  signal  with  DOA  0.  Equation  (1)  represents  the 

basic  form  of  the  theoretical  (uncompensated)  steering  vector. 

exp(-jk(Xj  cos  6  +  y,  sin  6)) 
exp(-jk(x2  cos  0  -H  y2  sin  0)) 
exp(-jk(x3  cos  0  +  y3  sin  0)) 

•  •  • 

exp(-jk(xM  COS0-I-  y^  sin0)) 


au(0)  = 


(1) 


In  the  above  equation,  k  is  the  wavenumber  of  the  incident  electromagnetic  radiation  that 
can  be  expressed  as  27t/A,.  An  array  manifold  matrix,  A  contains  as  its  columns  a 
collection  of  steering  vectors  corresponding  to  a  set  of  angles. 

DOA  algorithms  directly  use  knowledge  of  the  array  steering  vector  to  generate  a 
spatial  spectrum  that  gives  signal  power  versus  direction.  A  simple  approach  to  this 
DOA  estimation  is  to  assume  the  ideal  steering  vector  given  in  Equation  (1).  This  ideal 
steering  vector  is  a  function  of  array  geometry  and  incident  angle,  but  does  not  take  into 
account  mutual  coupling  effects.  Specifically,  Equation  (1)  does  not  take  into  account 
the  retransmission  of  energy  from  each  antenna  element  to  the  others. 

The  concept  of  a  spatial  signature  will  be  relevant  for  the  discussion  of  downlink 
beamforming  algorithms  and  is  introduced  here.  The  spatial  signature  model  considers 
uplink  antenna  array  output  due  to  a  single  user  to  be  a  linear  combination  of  the  steering 
vectors  of  all  the  direct  path  and  multipath  components. 

x(t)  =  a(0,  >,  (t)  +  ^aia(ei  )s,  (t)  =  v,s,  (t)  (2) 

i=2 

V,  is  defined  to  be  the  spatial  signature  of  mobile  user  1.  a(0,)s,(t)  is  the  direct  path 
component  and  aia(0j)s,(t) ,  2  <  i  <  Nm,  are  the  multipath  components.  The  oCj  contain 
the  relative  amplitude  and  phase  of  the  i*  signal  component  (modeling  relative 
attenuation  and  delay  respectively). 

B.  Distortion  Matrix 

A  distortion  matrix,  C,  is  used  to  encapsulate  the  effect  of  mutual  coupling  as  well 
as  the  amplitude  and  phase  distortions  caused  by  imperfect  antenna  array  elements.  This 


matrix,  which  can  be  estimated  from  collected  measurement  data  [1-5],  is  applied  to  the 
equation  of  the  theoretical  steering  vector  form  of  (1)  to  develop  a  steering  vector  that 
takes  mutual  coupling  and  unknown  sensor  gains  and  phases  into  aecount. 

ac(0)  =  C(0)a„(0)  (3) 

In  the  literature,  to  simplify  the  problem,  the  distortion  matrix  is  generally 
considered  to  be  independent  of  angle.  This  assumption  will  be  considered  in  Section  IV . 
Array  calibration  methods  attempt  to  estimate  the  matrix  C  algorithmically  off-line  [2,3], 
on-line  [4,12],  or  speeifically  measure  ac(0)  using  field  measurements.  However,  one  of 
the  main  problems  in  estimating  the  distortion  matrix  is  that  it  is  hard  to  separate 
coupling,  gain,  and  phase  issues  from  environmental  faetors  such  as  tower  platform 
effects  and  other  scatterers  located  close  to  the  array  [9]. 

C.  Uplink  Direction  of  Arrival  Algorithms 

An  important  quantity  in  many  array  signal  processing  algorithms  is  the  Spatial 
Covariance  Matrix,  R.  It  is  defined  in  Equation  (4). 

R  =  E(x(t)x«(t))  (4) 

x”(t)  indicates  the  complex  conjugate  transpose  of  vector  x(t).  In  practice,  this  matrix 
is  estimated  using  N  snapshots  of  the  actual  antenna  array  output  (with  sampling  interval 
T),  as  shown  in  Equation  (5). 


(5) 


An  eigenvalue  decomposition  of  R  can  be  used  to  find  the  orthogonal  projector 
onto  the  estimated  noise  subspace,  n"*".  Using  this  information,  the  equations  for  the 
Bartlett  and  MUSIC  spatial  spectra  can  be  derived  [13]: 


BARTLETT 


a"(e)Ra(e) 

a“(6We) 


“  MUSIC 


^  ^  a“(e)n^a(e) 


(6) 

(7) 


D.  Downlink  Beamforming  Methods 

DOA  information  gathered  during  uplink  is  used  by  a  given  downlink 
beamforming  algorithm  for  the  purpose  of  transmitting  information  most  efficiently  from 
the  base  station  to  the  mobile  user.  Each  downlink  beamforming  method  has  its  own 
technique  for  how  best  to  accomplish  this  efficient  transmission  [14,15]. 

•  Dominant  DOA  (DomDOA)  method  -  For  a  given  mobile  user,  this  method 
takes  the  angle  at  which  there  is  the  greatest  received  power  during  uplink 
(determined  by  uplink  spatial  spectrum)  and  focuses  all  transmitted  energy  in 
that  direction. 

•  Pseudoinverse  DOA  (PseDOA)  method  -  In  addition  to  transmitting  in  the 
direction  of  the  strongest  DOA  of  the  desired  user  like  the  DomDOA  method, 
nulls  are  placed  in  the  antenna  radiation  pattern  at  both  the  non-dominant 
DOAs  of  the  desired  user  and  all  the  DOAs  of  any  other  mobile  users  in  the 
system. 

•  Spatial  Signature  (SS)  method  -  This  method  uses  the  spatial  signature  of  the 
mobile  user  when  constructing  the  downlink  weight  vector.  It  differs  from  the 


prior  two  methods  in  that  instead  of  focusing  signal  energy  only  in  the 
direction  of  the  dominant  DO  A  of  the  desired  mobile  user,  it  steers  the  beam 
of  signal  energy  in  all  directions  corresponding  to  the  desired  mobile  user. 

•  Pseudoinverse  Spatial  Signature  (PseSS)  method  -  Whereas  the  SS  method 
tried  to  maximize  SNR  without  inserting  nulls,  this  method  maximizes  SNR 
while  inserting  nulls.  Specifically,  nulls  are  inserted  at  the  DO  As 
corresponding  to  the  spatial  signature  of  the  other  mobile  users  in  the  system. 

It  is  important  to  note  that  the  SS  and  PseSS  methods  do  not  use  steering  vectors 
at  all  when  constructing  downlink  beamforming  weights.  While  the  spatial  signature  can 
be  modeled  as  a  linear  combination  of  steering  vectors,  this  decomposition  is  unnecessary 
when  determining  the  spatial  signature  and  using  spatial  signature  based  downlink 
beamforming  algorithms. 

III.  Mutual  Coupling  Compensation 

Our  approach  to  mutual  coupling  compensation  entails  developing  an  estimate  for 
the  actual  array  response  5^(0)  given  an  arbitrary  geometry  antenna  array  via  MoM 

code.  The  MoM  code  used  in  this  study  is  NEC  (Numerical  Electromagnetics  Code), 
which  was  originally  developed  at  the  Lawrence  Livermore  Laboratory  [16]  to  perform 
MoM  analysis  to  model  the  interaction  between  electromagnetic  fields  and  wire 
segments.  A  wire  segment  model  of  an  arbitrary  geometry  antenna  array  can  be  specified 
in  an  input  file  to  NEC.  NEC  results  are  well  accepted  in  the  literature  and  versions  of 
NEC  are  freely  distributed  on  the  World  Wide  Web. 


Fundamentally,  the  MoM  represents  the  unknown  induced  current  on  an  object  in 
terms  of  a  given  set  of  basis  functions  and  enforces  Maxwell’s  boundary  conditions  at  a 
finite  number  of  points  on  the  object  being  modeled  [11].  Thus,  in  the  MoM,  each 
antenna  array  element  is  represented  as  a  wire  of  finite  thickness  divided  up  into 
segments.  At  the  heart  of  the  MoM  is  the  calculation  of  a  system  impedance  matrix  Zqq 
giving  the  coupling  between  segment  i  and  j  in  the  antenna  array  model  (1<  i,j  ^  Q) 
where  Q  is  the  total  number  of  wire  segments  in  the  model.  The  value  of  Q  needs  to  be 
chosen  to  be  large  enough  so  that  the  obtained  current  solution  on  each  of  the  antenna 
elements  converge  to  a  fixed  function.  In  general,  the  higher  the  number  of  segments,  the 
greater  the  validity  and  resolution  of  the  answer  (at  the  expense  of  greater  computational 
effort).  A  detailed  formulation  of  the  modeling  of  an  antenna  array  using  the  MoM  can 
be  found  in  the  work  of  Pasala  [7]  and  Adve  [10]. 

A  7-element  circular  array  was  used  in  simulations  and  measurements  in  this 
study.  The  simulated  antenna  elements  are  each  coaxial  dipole  elements  containing  four 
collinear  U2  dipole  antennas  (illustrated  in  Figure  1)  with  an  operating  frequency  of  1.88 
GHz.  Each  of  the  dipole  elements  is  represented  as  a  wire  divided  up  into  segments  for 
MoM  analysis.  Lump  loading  is  used  to  model  the  isolation  between  the  dipole  antennas 
and  the  load  impedance  of  the  array  elements. 

Two  types  of  mutual  coupling  compensation  using  MoM  were  considered  in  this 
study.  In  the  first  type  of  compensation,  MoM  calculations  were  used  to  determine  the 
true  steering  vectors  in  the  azimuth  plane  perpendicular  to  the  axis  of  the  array  elements. 
This  complete  array  manifold  is  then  used  when  performing  uplink  DOA  analysis.  In  the 


second  compensation  technique,  based  upon  [1-5,9],  the  steering  vectors  at  only  a  few 
angles  are  computed,  and  Equation  (3)  is  solved  for  the  distortion  matrix  (assuming 
angular  independence)  using  the  equation: 

c  =  a„a1(a,„aL)"'  (8) 

where  the  columns  of  Atme  and  Atheo  are  the  MoM  determined  (compensated)  steering 
vectors  and  the  theoretical  (uncompensated)  steering  vectors  respectively,  corresponding 
to  a  particular  set  of  angles.  This  distortion  matrix  is  then  applied  using  Equation  (3)  to 
determine  compensated  steering  vectors.  Field  measurement  comparsion  of  the  results  of 
these  two  types  of  compensation  allow  the  assumption  of  angular  independence  of  the 
distortion  matrix  to  be  examined. 


IV.  Results 

A.  Steering  Vector  Relative  Angle  Change 

We  first  consider  the  difference  between  the  ideal  (from  Equation  (1))  and 
compensated  (from  NEC  simulations)  array  steering  vectors.  As  a  metric,  we  use  the 
calculation  of  relative  angle  change  (RAC)  between  theoretical  and  compensated  array 
vectors.  For  a  compensated  steering  vector,  ac,  and  an  uncompensated  steering  vector, 
au,  this  metric  can  be  written  as: 


RAC(%)  =  lOOx 


(9) 


This  metric  has  been  used  in  simulations  and  measurements  [15,17]  to  provide  a 
measure  of  the  difference  between  spatial  signatures.  It  is  used  in  this  study  to  quantify 
the  difference  between  compensated  and  uncompensated  steering  vectors,  and  is  shown 


in  Figure  2.  This  graph  shows  that  there  is  an  appreciable  difference  between  the 
compensated  and  uncompensated  steering  vectors. 

B.  Experimental  Setup 

Data  collection  was  performed  outdoors  on  an  open  field  in  the  Pickle  Research 
Campus  at  the  University  of  Texas  at  Austin.  There  were  two  signal  generators  used  to 
represent  mobile  users  transmitting  to  the  smart  antenna  basestation.  The  first  signal 
generator  was  held  stationary.  The  second  signal  generator  was  moved  to  five  distinct 
angles  at  approximately  the  same  distance  and  the  same  transmission  power  as  the  first 
signal  generator.  Data  was  then  collected  using  the  smart  antenna  testbed  and  uplink 
spatial  spectra  were  considered  for  each  generator  transmitting  individually  and  both 
generators  transmitting  together.  The  position  of  each  of  the  signal  generator  locations  is 
listed  in  Table  1. 

C.  Uplink  Direction  of  Arrival  Algorithm 

Figure  3  shows  the  spatial  spectra  generated  using  uncompensated  and 
compensated  steering  vectors  for  measurement  data  due  to  a  single  mobile  user  (referred 
to  as  “Stationary  User”  in  Table  1).  In  this  figure,  we  consider  the  application  of  the  full 
array  manifold  calculated  using  the  MoM.  Figure  3a  contains  Bartlett  DO  A  algorithm 
output  and  Figure  3b  contains  the  MUSIC  spatial  spectrum.  These  spatial  spectra  show 
that  while  the  effect  on  the  detected  DOA  is  minimal  for  both  uncompensated  and 
compensated  results,  sidelobe  levels  are  significantly  reduced  through  use  of  mutual 
coupling  compensation.  Specifically,  Figure  3b  shows  the  t5^ical  (for  all  tested  cases) 


sidelobe  reduction  of  approximately  3  dB.  While  compensation  has  little  effect  on  the 
main  lobe  of  the  Bartlett  spatial  spectra  in  Figure  3  a,  the  effect  is  more  pronounced  on  the 
MUSIC  spatial  spectra  in  Figure  3b. 

Figure  4  shows  the  same  results  as  Figure  3  for  a  different  mobile  user. 
Specifically,  it  represents  the  Bartlett  and  MUSIC  DOA  algorithm  output  given  the  full 
NEC  calculated  array  manifold  for  a  single  mobile  user  (referred  to  as  “Position  1”  in 
Table  1).  Again,  this  figure  shows  that  mutual  coupling  compensation  has  a  greater 
effect  on  the  main  lobe  of  MUSIC  spatial  spectrum  output  shown  in  Figure  4b  compared 
with  the  effect  on  the  main  lobe  of  the  Bartlett  spatial  spectrum  shown  in  Figure  4a.  Both 
spatial  spectra  in  Figure  4  show  a  reduction  in  unwanted  sidelobe  levels. 

Figure  5  shows  the  uncompensated  and  compensated  MUSIC  spatial  spectra  for  a 
situation  in  which  there  was  two  mobile  users  (“Stationary  User”  and  “Position  5”  in 
Table  1)  transmitting  simultaneously.  This  figure  again  illustrates  the  previous  results  of 
significant  sidelobe  reduction  along  with  sharpening  of  the  lobes  corresponding  to  mobile 
users.  Specifically,  note  that  the  DOA  of  the  stationary  user  is  made  much  more 
prominent  compared  to  the  sidelobes  due  to  compensation.  These  sidelobes  could  be 
mistaken  for  either  another  mobile  user  or  multipath  signal  energy  -  either  of  which  has 
adverse  consequences  for  downlink  beamforming  methods  such  as  the  PseDOA  method. 
In  general  for  all  tested  cases  (single  and  multiple  user  scenarios),  mutual  coupling 
compensation  increased  lobes  in  the  spatial  spectra  corresponding  to  desired  users’  DO  As 
and  decreased  all  other  sidelobes. 

Figure  6  shows  the  effect  of  assuming  the  distortion  matrix  is  independent  of 
angle  using  the  same  spatial  spectrum  considered  in  Figure  3b.  For  all  of  the  results  in 


the  paper  not  including  those  of  Figure  6,  the  compensation  technique  involved  using  the 
entire  array  manifold  computed  through  NEC.  In  Figure  6,  this  technique  is  referred  to  as 
the  “Look-up  Table”  approach.  The  other  tested  compensation  technique  is  to  compute 
the  distortion  matrix  (assuming  angular  independence)  using  Equation  (8)  with  a 
relatively  small  number  of  MoM-generated  and  theoretical  steering  vectors.  Results  from 
this  study  show  that  making  use  of  the  assumption  that  C  is  independent  of  angle  yields  a 
performance  increase  over  theoretical  steering  vector  use.  This  performance  increase 
does  not  quite  match  the  benefit  due  to  the  use  of  the  complete  MoM-generated  array 
manifold.  However,  it  should  be  pointed  out  that  using  the  “Look-up  Table”  approach 
has  the  disadvantage  of  requiring  more  storage  than  a  computed  distortion  matrix. 

D.  Downlink  Beamforming  Algorithm 

To  determine  the  effect  that  mutual  coupling  compensation  has  on  downlink 
beamforming,  the  MoM  model  for  the  antenna  array  was  used.  The  uplink  spatial  spectra 
from  measurement  data  for  all  of  the  tested  cases  were  used  to  compute  downlink 
beamforming  weights  for  scenarios  involving  the  Stationary  User  (User  1)  and  the  mobile 
user  at  each  of  5  different  tested  positions  (User  2).  Figure  7  shows  the  resulting 
simulated  radiation  pattern  from  using  each  of  the  four  downlink  beamforming  methods 
described  in  Section  II-D.  Figure  7a  shows  the  downlink  radiation  pattern  for  the 
Stationary  User  and  Figure  7b  shows  the  pattern  for  the  mobile  user  at  Position  4.  This 
figure  shows,  for  the  DOA  based  methods  DomDOA  and  PseDOA,  the  effects  of  using 
compensated  and  uncompensated  steering  vectors.  As  seen  in  the  top  pattern  of  Figures 
7a  and  7b,  use  of  compensated  steering  vectors  does  not  significantly  improve  the 


performance  of  the  DomDOA  method.  Specifically,  both  of  these  patterns  show  that  the 
main  lobes  of  the  uncompensated  and  compensated  DomDOA  method  are  practically  the 
same.  Furthermore,  both  patterns  show  similar  performance  to  that  of  the  SS  method. 
However,  the  bottom  pattern  of  Figures  7a  and  7b  shows  that  there  is  significant 
improvement  in  the  performance  of  PseDOA  method  when  using  compensated  steering 
vectors.  These  patterns  clearly  show  an  improved  main  lobe  when  using  the 
compensated  PseDOA  method  instead  of  the  uncompensated  PseDOA  method.  In 
addition,  these  patterns  show  that  use  of  compensation  allows  the  PseDOA  method  to 
work  comparably  to  the  PseSS  method. 

To  quantify  the  results  illustrated  by  Figure  7,  Table  2  shows  the  gain  in  signal 
power  for  each  mobile  user  that  is  realized  through  use  of  the  compensated  PseDOA 
method  over  the  uncompensated  PseDOA  method  and  the  PseSS  method.  Signal  power 
is  determined  by  the  value  of  the  radiation  plot  at  the  known  DO  A  of  the  desired  user. 
As  seen  in  this  table,  the  use  of  compensated  steering  vectors  instead  of  uncompensated 
vectors  in  the  PseDOA  method  (Column  1  for  User  1  and  Column  3  for  User  2)  nearly 
always  significantly  improves  transmitted  signal  power.  Furthermore,  this  table  shows 
that  use  of  compensated  steering  vectors  allows  the  PseDOA  method  to  work  just  as  well 
as  the  PseSS  method  (Column  2  for  User  1  and  Column  4  for  User  2),  which  is  normally 
not  true  [14]. 

Table  3  shows  the  improvement  in  SIR  for  the  two  mobile  users  using  the 
compensated  PseDOA  method  over  the  uncompensated  PseDOA  method.  In  this  table, 
signal  power  is  determined  as  before  and  interference  power  for  one  mobile  user  is 
determined  by  the  value  at  that  user’s  DO  A  on  the  other  user’s  radiation  plot.  This  table 


shows  that  use  of  compensated  steering  vectors  instead  of  uncompensated  vectors  often 
significantly  improves  SIR  (Column  1  versus  Column  2  for  User  1  and  Column  3  versus 
Column  4  for  User  2).  The  results  in  this  table  also  show  that  SIR  performance  is  more 
stable  as  a  function  of  user  position  when  using  compensated  steering  vectors  (Column  1 
for  User  1  and  Column  3  for  User  2).  This  result  has  beneficial  implications  for  the 
minimization  of  cellular  call  dropping. 

V.  Conclusion 

In  this  paper,  we  have  examined  the  benefit  of  mutual  coupling  compensation 
during  uplink  DOA  estimation  and  during  downlink  beamforming.  Results  from  our 
measurement  campaign  allow  us  to  make  several  comments  regarding  mutual  coupling 
compensation.  It  has  been  shown  through  the  use  of  MoM  calculations  that  steering 
vectors  taking  into  account  mutual  coupling  differ  significantly,  in  terms  of  relative  angle 
change,  from  theoretical  steering  vectors.  By  using  more  accurate  steering  vectors,  we 
observed  the  expected  performance  improvement.  The  performance  improvement  in 
terms  of  direction  finding  comes  from  reduced  sidelobe  levels,  which  allows  for  easier 
distinction  of  signal  sources.  The  DOA  performance  increase  also  translates  into 
improved  downlink  SIR  performance,  allowing  a  DOA-based  downlink  beamforming 
method  (Pseudoinverse  DOA  method)  to  work  comparably  to  spatial  signature  based 
methods.  This  work  showed  how  off-line  calculations  using  basic  electromagnetic 
computations  could  lead  to  significant  improvements  in  smart  antenna  system 
performance. 
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Figure  3  -  Experimental  Data  Stationary  User  Spatial  Spectra  -  (A)  Bartlett  DOA 
Algorithm  (B)  MUSIC  DOA  Algorithm 
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Mobile  Position  5 
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ABSTRACT—  An  approximate-inverse  preconditioner  based  on  the  pre-defined  wavelet 
packet  (PWP)  basis  is  proposed  for  the  fast  iterative  solution  of  electromagnetic  integral 
equations.  The  PWP  basis  is  designed  to  achieve  a  sparse  representation  of  the  moment 
matrix  and  the  preconditioner  is  constructed  by  inverting  the  block-diagonal 
approximation  of  the  PWP-based  moment  matrix  and  transfomung  the  results  into  the 
space  domain.  Numerical  results  show  that  the  PWP  preconditioner  is  effective  in 
accelerating  the  convergence  rate  of  iterative  solution  to  moment  equations.  It  is  also 
demonstrated  that  by  properly  designing  the  block-diagonal  matrix  and  computing  the 
matrix  elements,  the  total  computational  complexity  and  memory  costs  for  the 
preconditioner  can  be  kept  to  OG'WogN). 

Index  Terms  —  Computational  Electromagnetics,  Preconditioning,  Wavelet  Packet  Basis. 

1.  Introduction 

Iterative  algorithms  are  commonly  used  to  solve  large-scale  moment  equations 
resulting  from  electromagnetic  integral  equations.  The  computational  cost  of  iterative 
solution  is  proportional  to  both  the  moment  matrix-vector  multiplication  operation  and 


the  number  of  iterations  required  for  a  convergent  solution.  The  multilevel  fast  multipole 
method  (MLFMM)  has  been  developed  to  reduce  the  time  complexity  and  memory  of  the 
dense  matrix-vector  multiplication  operation  from  O(N^)  to  O(NlogN)  [1-3].  However, 
the  total  solution  time  remains  dependent  on  the  number  of  iterations  required  to  achieve 
an  accurate  solution.  If  the  scatterer  contains  near-resonant  structures  such  as  deep 
cavities,  the  moment  matrix  is  not  well  conditioned  and  the  convergence  rate  of  the 
iterative  solution  can  become  very  slow  [4,  5].  In  these  cases,  a  preconditioning  matrix 
[P]  is  desired  to  accelerate  the  convergence  rate.  We  consider  a  left-preconditioning 
system: 

[P][Z]J  =  [P]E  (1) 

where  [Z]  is  the  moment  matrix,  J  is  the  unknown  induced  current  vector  and  E  is  the 
excitation  vector.  To  achieve  the  best  preconditioning,  [P]  should  be  as  close  to  the 
inverse  of  [Z]  as  possible.  Hence,  [P]  is  called  an  approximate-inverse  preconditioner.  In 
this  work,  we  set  out  to  construct  an  approximate-inverse  preconditioner  [P]  that  satisfies 
the  following  conditions: 

(a)  [P]  is  an  accurate  approximation  to  [Z]'*  to  be  an  effective  preconditioner. 

(b)  The  total  computational  complexity  and  memory  requirement  to  construct  the 
preconditioner  are  of  O(NlogN). 

(c)  The  total  computational  complexity  and  memory  requirement  to  implement  the 
preconditioning  operation  in  (1)  are  of  O(NlogN). 

Conditions  (b)  and  (c)  are  set  to  keep  the  complexity  of  the  preconditioning  on  par  with 
that  of  the  MLFMM,  making  our  preconditioner  applicable  for  the  MLFMM.  The 
traditional  preconditioning  methods,  such  as  incomplete  LU  factorization  (ILU), 


Frobenius-norm-based  approximate  inverse  and  polynomial  preconditioners  could  be 
effective  [6,  7].  But  they  are  computationally  too  complicated  to  meet  conditions  (b) 
and  (c).  In  [5],  an  approximate-inverse  preconditioner  was  constructed  by  inverting  a 
block-banded  form  of  the  original  moment  matrix.  Low  complexity  was  maintained  by 
choosing  a  bandwidth  that  is  independent  of  problem  size. 

Our  approach  to  the  problem  is  to  derive  an  effective  preconditioner  meeting 
conditions  (a)-(c)  based  on  the  wavelet  packet  basis.  With  the  conventional  subsectional 
basis,  [Z]  in  (1)  is  dense,  making  it  difficult  to  find  an  effective  approximate-inverse 
preconditioner.  However,  if  we  can  transform  the  moment  equation  to  a  new  basis  to 
make  the  transformed  moment  matrix  sparse  and  diagonal  dominant,  it  will  be  much 
easier  to  find  an  effective  preconditioner.  We  have  recently  proposed  the  Pre-defined 
Wavelet  Packet  (PWP)  basis  for  the  efficient  representation  of  moment  matrices  [8,  9]. 
The  PWP  basis  was  designed  to  match  the  oscillatory  nature  of  the  Green’s  kernel  in  the 
integral  equation.  Numerical  results  showed  that  the  PWP-transformed  moment  matrices 
are  strongly  diagonal  dominant,  and  have  about  C>(NlogN)  non-zero  elements  for  large- 
scale  problems.  In  [10],  a  preconditioner  for  the  PWP-transformed  moment  matrix  was 
constructed  to  accelerate  the  convergence  of  the  transformed  matrix  equation.  In  this 
paper,  we  shall  demonstrate  that  an  effective  approximate-inverse  preconditioner  in  the 
space  domain  can  be  derived  from  the  PWP  basis  with  no  more  than  O(NlogN) 
computational  or  memory  cost.  It  is  worthwhile  to  point  out  that  there  have  been  some 
recent  works  on  using  the  conventional  wavelet  transform  to  precondition  linear  systems 
for  iterative  solutions  [11-14].  However,  the  requirements  for  computational  cost  and 
memory  are  at  least  (9(N^)  to  implement  the  wavelet  transform  and  the  preconditioning 


operation.  In  this  work,  we  overcome  the  6>(N^)  complexity  and  memory  bottlenecks  by 
computing  the  PWP  preconditioner  directly  from  the  PWP  basis  functions  [9]. 

This  paper  is  organized  as  follows.  In  Section  H,  we  give  a  brief  review  of  the 
PWP  basis  and  its  representation  of  the  moment  matrices.  In  Section  HI,  we  describe  the 
design  and  construction  of  our  approximate-inverse  preconditioner  in  the  space  domain 
from  the  sparse  PWP-based  moment  matrix.  We  show  that  the  total  cost  for  the 
construction  and  preconditioning  operation  can  be  kept  under  O(NlogN).  Numerical 
results  are  presented  in  Section  FV.  Finally  we  draw  some  conclusions  in  Section  V. 


11.  PWP  Basis  and  Its  Representation  of  Moment  Equations 

Wavelet  packet  basis  is  a  set  of  orthonormal  functions  derived  from  the  shift  and 
dilation  of  a  basic  wavelet  function,  and  can  be  considered  as  a  generalization  of  the 
conventional  wavelet  basis  [15].  From  the  point  of  view  of  function  space 
decomposition,  the  wavelet  packet  basis  space  is  generated  from  the  decomposition  of 
both  the  scaling  function  space  and  the  corresponding  wavelet  function  space.  Let  us 
assume  that  \i/(x)  is  the  wavelet  packet  basis  function  with  the  finest  spatial  resolution 
available  for  signal  representation.  Using  the  “two-scale  equations,”  we  can  express  the 
wavelet  packet  basis  functions  in  the  next  scale  as  [16]: 


yfl{n)  =  Y,yf{k)h{2n-k) 
k 

W\  (.n)  =  X  yf(k)g(2n  -  k) 

it 


(2) 


where  {h(k)}  and  {g(k)}  are  the  impulse  responses  of  two  quadrature  filters  H  (low-pass) 
and  G  (high-pass),  respectively.  The  functions  in  the  next  scale  become  coarser  in  spatial 
resolution  and  finer  in  spectral  resolution  after  the  filtering  and  down-sampling  in  (2). 


The  same  procedure  can  be  applied  recursively  to  the  outputs  of  (2)  into  subsequent 
scales.  Conversely,  the  decomposition  results  in  (2)  can  also  be  used  to  reconstruct  the 
original  sequence  by  using  another  pair  of  quadrature  filters  P  and  Q,  which  are  the 
mirror  filters  of  H  and  G,  respectively.  This  reconstruction  procedure  can  be  expressed 
as: 

y/ix)  =  X  (*)  +  S  9  (^  -  (*)  (3) 

k  k 

where  {p(k)}  and  {q(k)}  are  the  impulse  responses  of  P  (low-pass)  and  Q  (high-pass), 
respectively.  The  functions  become  finer  in  spatial  resolution  and  coarser  in  spectral 
resolution  through  filtering  and  up-sampling  in  (3). 

A  complete  and  orthogonal  wavelet  packet  basis  can  be  generated  from  a  spectral 
decomposition  tree  that  starts  from  an  initial  mother  function  \|/(x)  with  the  finest  spatial 
resolution  by  using  recursive  two-channel  filtering  and  down-sampling  in  (2).  It  can  be 
shown  that  the  decomposed  functions  at  the  outermost  branches  of  the  tree  satisfy 
orthogonality  and  completeness  for  any  decomposition  trees,  and  thus  constitute  a 
wavelet  packet  basis  set  [17].  In  this  context,  the  conventional  wavelet  transform  (CWT) 
basis  can  be  considered  as  a  special  case  of  the  wavelet  packet  basis  when  the 
decomposition  tree  grows  along  only  the  lowest  frequency  branch.  However,  we  have 
shown  in  [18]  that  the  structure  of  the  CWT  tree  is  not  optimal  for  representing  the 
moment  matrix  from  electrodynamic  integral  equations.  Instead,  we  have  proposed  a 
special  wavelet  packet  tree  for  electrodynamic  problems  [8,  9,].  This  class  of  wavelet 
packet  bases,  which  we  term  the  pre-defined  wavelet  packet  (PWP)  basis,  has  a  spectral 
decomposition  tree  that  zooms  in  along  the  free-space  wave  number  ko.  The  motivation 
for  this  tree  structure  is  to  ensure  that  the  basis  is  well-matched  to  the  oscillatory  nature 


of  the  Green’s  function  kernel  in  the  integral  equation.  Figs.  1(a)  and  1(b)  show 
respectively  the  conventional  wavelet  decomposition  tree  and  the  PWP  decomposition 
tree  for  N=32.  In  the  PWP  tree,  the  center  frequency  of  its  deepest  branch  is  as  close  as 
possible  to  ko.  Detailed  discussion  on  the  construction  of  the  PWP  basis  can  be  found  in 
[9]. 

Once  the  PWP  basis  is  defined,  it  is  straightforward  to  transform  the  original 
moment  equation  into  the  new  basis  representation  as  follows: 


[Z]J  =  E 

(4) 

where 

[Z]  =  [M]^[Z][M] 

(5) 

J  =  [Mfj 

(6) 

E=[MfE 

(7) 

[M]^  denotes  the  transpose  of  [M],  and  is  the  PWP  transformation  matrix.  It  changes  the 
original  expansion/testing  functions  from  the  standard  subsectional  bases  into  the  PWP 
bases.  Note  that  [M]^  is  equal  to  [M]"'  since  [M]  is  an  orthonormal  transformation. 

[Z],  J,  and  E  are  the  moment  matrix,  induced  current,  and  excitation  vector 
represented  in  the  PWP  basis,  respectively.  Our  numerical  results  from  [9]  showed  that 
the  transformed  moment  matrix  under  the  PWP  basis  is  very  sparse.  In  particular  the 
number  of  above-threshold  elements  in  the  transformed  moment  matrix  grows  at  a  rate  of 
(9(n'’^)  for  small-sized  problems,  and  approaches  (9(NlogN)  for  large-scale  problems. 

Furthermore,  the  significant  elements  in  [Z]  are  located  mainly  along  the  diagonal  or 
near-diagonal  positions.  Unfortunately,  the  exact  locations  of  the  significant  elements  are 
difficult  to  predict.  Consequently,  the  complexity  of  the  algorithm  in  solving  the 
complete  moment  equation  is  still  hampered  by  an  computational  bottleneck,  since 


every  element  of  the  matrix  must  be  computed  to  determine  whether  it  is  small  enough  to 
be  discarded.  However,  for  the  proposed  preconditioning  application  in  this  paper,  we 
shall  show  in  the  next  section  that  an  effective  approximate-inverse  preconditioner  can  be 
constructed  by  computing  only  those  elements  within  the  diagonal  blocks. 

III.  PWP-Based  Approximate-Inverse  Preconditioner  in  the  Space  Domain 

We  will  first  outline  our  approach  to  constructing  the  PWP-based  preconditioner. 
Since  the  PWP  transformed  moment  matrix  [Z]  is  very  sparse  and  diagonal  dominant, 
we  approximate  [  Z  ]  by  a  block-diagonal  matrix  [  Z  bd]  that  consists  of  the  near-diagonal 
terms  of  [Z].  The  approximate-inverse  preconditioner  [P]  defined  in  (1)  is  then 

constructed  by  inverting  the  block-diagonal  matrix  [  Z  bd]  and  transforming  the  resulting 
matrix  back  to  the  space  domain: 

[P]  =  [M][Zbdr'[Mf  (8) 

Therefore,  the  preconditioned  moment  equation  becomes: 

[M]  [  Z  bd]  ‘  [M]^[Z]  J  =  [M]  [  Z  bd]'‘  [M]^E  (9) 

The  inverse  of  the  block  diagonal  matrix  [Zbdl  is  simply  the  inverse  of  its  individual 
diagonal  blocks,  and  the  inverted  matrix  [Zbdl  *  remains  block  diagonal.  By  properly 
choosing  the  block  sizes  of  [Zbdl.  we  can  limit  the  computational  cost  of  the  inversion 
while  maintaining  the  effectiveness  of  the  preconditioner. 

The  actual  implementation  of  the  preconditioner  will  now  be  considered.  If  we 
solve  equation  (9)  using  an  iterative  algorithm  such  as  the  conjugate  gradient  method,  the 
complexity  of  each  iteration  is  proportional  to  that  of  the  product  between  the  matrix 
[M],  [Zbdr‘,  [M]  ^  or  [Z]  and  a  vector.  The  product  [Z]  and  a  vector  can  be  implemented 


using  the  MLFMM  algorithm  with  O(NlogN)  complexity  and  memory  requirement. 
Because  the  PWP  transformation  matrix  [M]^  has  about  O(NlogN)  non-zero  elements, 
the  multiplication  of  [M]^  or  [M]  with  a  vector  is  also  of  C>(NlogN)  in  complexity. 
Therefore,  if  we  can  limit  both  the  number  of  non-zero  elements  in  [  Z  bd] '  and  the 
complexity  of  constructing  [Zbd]  *  to  O(NlogN),  the  total  cost  per  iteration  in  (9)  will 
remain  at  (9(NlogN). 

Our  design  of  the  block-diagonal  matrix  [  Z  bd]  is  shown  in  Fig.  2.  The  blocks  are 
designed  so  as  to  capture  the  most  significant  elements  in  the  PWP-transformed  moment 
matrix  [Z].  The  block-diagonal  matrix  consists  of  the  following:  (i)  a  block  centered 
about  the  spectral  frequency  ko  with  block  size  Mi,  (ii)  a  set  of  diagonal  blocks  in  the 
remaining  upper-left  region  with  block  size  M2,  and  (iii)  a  set  of  diagonal  blocks  in  the 
lower-right  region  with  block  size  M3.  The  number  of  diagonal  blocks  with  size  M2  is 
[(N/2-  Ml)/  M2],  and  that  with  size  M3  is  [(N/2)/  M3].  The  block  from  (i)  is  usually  the 
densest  part  of  [  Z  ],  since  its  elements  correspond  to  the  interactions  between  PWP  bases 
with  the  longest  spatial  extent  and  spectral  frequency  close  to  ko.  They  tend  to  have  the 
strongest  interaction  with  each  other.  The  blocks  from  (iii),  on  the  other  hand,  are  nearly 
diagonal,  since  their  elements  correspond  to  the  interactions  between  PWP  bases  with  the 
highest  spectral  frequency  and  the  shortest  spatial  extent.  Therefore,  we  generally  choose 
block  sizes  with  Mi  >  M2>  M3.  Furthermore,  as  N  is  increased  we  keep  the  sizes  M2  and 
M3  fixed,  but  let  Mi  grow  slightly  with  problem  size  at  a  rate  of  to  capture  even  more 
of  the  dominant  elements  in  [Z].  This  growth  rate  is  chosen  to  ensure  that  the  cost  of 
inverting  the  block  remains  bounded  by  N,  as  the  cost  of  the  inversion  is  proportional  to 
the  cube  of  the  matrix  dimension.  The  inverse  of  [Zbd]  is  equivalent  to  the  inverse  of 


each  of  the  blocks  of  [Zbd]  in  Fig.  2.  Therefore,  under  our  construct,  the  total 
computational  complexity  to  invert  [  Z  bd]  is  proportional  to  N  and  the  resulting  number  of 
nonzero  elements  in  [ Zbd] '  is  also  proportional  to  N. 

Next,  we  discuss  the  computation  of  the  required  elements  in  [Zbd].  Each  element 
Z  (m,  n)  of  the  matrix  [  Z  ]  can  be  directly  expressed  as  [9]: 

Z(m,n)  =  YLy'n.(i)G{i,j)¥nU)  (10) 

i  j 

where  G(i,  j)  is  the  free  space  Green’s  function,  i  and  j  represent  the  indices  of  the 
discretization  grid  in  the  space  domain,  and  ij/k  is  the  k-th  PWP  basis  function.  It  can  be 
seen  from  (10)  that  the  computational  complexity  for  each  matrix  element  is  proportional 
to  the  product  of  the  spatial  supports  of  the  two  corresponding  basis  functions.  The 
support  of  a  PWP  basis  function  is  related  to  the  depth  of  the  branch  in  the  PWP 
decomposition  tree.  Let  us  assume  that  the  basis  functions  corresponding  to  branches  in 
the  first  stage  of  the  PWP  decomposition  tree  have  spatial  support  L  (i.e.,  the  order  of  the 
wavelet  filter).  The  complexity  for  computing  a  diagonal  (or  near-diagonal)  element  is 
then  I?.  The  support  of  a  PWP  basis  function  doubles  if  it  corresponds  to  the  branches 
on  the  next  deeper  stage  in  the  decomposition  tree  [9,  16].  Therefore,  the  spatial  support 
of  the  basis  functions  generated  from  the  tree  branches  at  stage  k  is  2'^'*L,  and  the 
complexity  to  compute  an  element  from  basis  functions  in  the  same  stage  is  (2  ‘  L) . 
Since  the  maximum  possible  depth  of  the  wavelet  decomposition  tree  is  Kmax=log2N,  the 
cost  for  computing  an  element  using  bases  from  this  maximum  depth  would  then  be  N  . 
However,  as  we  shall  discuss  in  the  next  section,  the  best  preconditioning  performance  is 
often  achieved  using  a  tree  with  depth  much  less  than  the  maximum  possible  depth. 
Similar  observation  has  also  been  reported  in  [11].  Therefore,  if  we  impose  that  the 


maximum  depth  index  of  the  PWP  decomposition  trees  is  a  small  number  K,  the  number 
of  operations  needed  to  compute  an  element  in  [Zbdl  will  be  between  L  and  4  L  . 
With  0(N)  elements  in  [Zbdl.  the  total  complexity  to  compute  the  matrix  [Zbdl  is  about 
C>(N)  operations. 

To  summarize,  we  have  designed  a  preconditioner  based  on  the  block-diagonal 
approximation  of  the  PWP-transformed  moment  matrix.  The  steps  to  construct  and  carry 
out  the  preconditioning  operation  entail:  (i)  computing  the  elements  in  the  diagonal 
blocks  of  [Zbdl  one  term  at  a  time  using  (10),  (ii)  inverting  [Zbdl  one  block  at  a  time  to 
arrive  at  its  inverse  [Zbdl’,  and  (iii)  carrying  out  the  preconditioning  operation  by  a 
successive  set  of  sparse  matrix-vector  multiplication  operations  in  (9).  It  is  shown  that 
the  computational  complexity  of  the  each  of  the  three  steps  can  be  limited  to  within 
C>(NlogN)  operations,  thus  making  the  proposed  preconditioner  compatible  with  the 
MLFMM  in  solving  moment  equations. 

IV.  Numerical  Results 

Two  two-dimensional  targets,  an  open-ended  inlet  and  a  bent  structure  shown  in 
Fig.  3,  are  chosen  to  test  our  PWP  preconditioner.  With  a  large  depth-to-opening  ratio 
(llm),  the  moment  matrices  from  the  scatterers  are  ill-conditioned,  and  the  convergence 
rate  is  slow  if  an  iterative  solver  is  used  to  solve  the  systems.  For  the  inlet  structure,  if  the 
thickness  parameter  t  is  equal  to  zero,  the  electric  field  integral  equation  (EFIE)  is  use  to 
construct  the  moment  equation.  Otherwise,  the  combined  field  integral  equation  (CFIE) 
is  used.  For  the  bent  structure,  EFIE  is  always  used  to  find  the  solution.  The 
discretization  of  the  structure  is  chosen  to  be  fixed  at  0.1  wavelength  while  the  electrical 


size  of  the  structure  is  varied.  In  all  cases  the  incident  plane  wave  is  E-polarized  and 
makes  an  angle  of  45°  with  respect  to  the  inlet  opening. 

Fig.  4  is  the  moment  matrix  [Z]  represented  by  the  PWP  basis  for  the  inlet 
structure  with  dimensions  15:1:1  and  N=256.  It  is  shown  in  logarithmic  scale 
without  any  thresholding.  We  find  that  the  most  significant  elements  in  the  transformed 
matrix  match  the  patterns  we  chose  for  the  PWP  preconditioner  in  Fig.  2.  The  densest 
portion  of  the  matrix  is  located  in  the  upper  left  region,  with  the  strongest  elements  near 
the  ko  spectral  position.  The  wavelet  filter  used  to  generate  the  PWP  basis  functions  is 
Daubechies  filter  with  order  16  (with  7  vanishing  moments). 

We  use  both  the  conjugate-gradient  squared  (CGS)  method  and  the  biconjugate 
gradient  stabilized  (BiCGSTAB)  method.  Both  of  them  are  effective  iterative  solvers  for 
linear  equations  with  nonsymmetric  coefficient  matrix  [19].  CGS  is  used  in  most  of  the 
cases  in  this  work  since  it  converges  faster  than  BiCGSTAB.  However,  we  use  the 
BiCGSTAB  to  compare  the  detailed  convergence  rates  of  the  different  preconditioners,  as 
the  convergence  behavior  of  the  CGS  is  not  very  smooth.  In  the  application  of  the 
iterative  solvers,  the  initial  current  guess  is  Jo=0,  and  the  stopping  criterion  used  is  when 
the  relative  residual  error  ||^‘„||/||^o||  - 10”*,  where  r„  =  E-[Z]Jn  and  ro=  E. 

We  first  examine  how  the  iteration  number  changes  as  the  electrical  size  of  the 
scatterer  is  increased  in  the  following  three  cases: 

Case  1:  The  scatterer  is  the  inlet  structure  shown  in  Fig.  3(a)  with  t  =  0  and 
/:m=15:l.  The  moment  equation  is  formulated  with  the  EFIE.  For  the  PWP 
preconditioner,  we  choose  M2=32,  M3=l,  and  Mi  grows  at  a  rate  of  \fN  starting  from  32 


atN=128. 


Case  2:  The  scatterer  is  the  inlet  structure  with  /:7n:t=  15:1:1.  The  moment 
equation  is  formulated  with  the  CFIE.  The  PWP  preconditioner  is  constructed  the  same 
way  as  that  in  Case  1 . 

Case  3:  The  scatterer  is  the  bent  structure  shown  in  Figure  3(b)  with  /:m=12:l. 
The  EFE  is  used,  and  the  parameters  for  the  PWP  preconditioner  are  chosen  as:  M2= 

M3=  32,  and  Mi  grows  at  a  rate  of  ifW  starting  from  32  at  N=128. 

The  iteration  numbers  required  versus  problem  sizes  are  shown  in  Figs.  5(a)-(c) 
for  Cases  1-3,  respectively.  We  increase  the  problem  size  by  proportionally  increasing 
the  electrical  size  of  the  scatterer  while  keeping  the  discretization  interval  fixed  at  0.1 
wavelength.  The  results  for  the  moment  equation  without  any  preconditioning  and  that 
with  the  PWP  preconditioner  are  plotted.  For  comparison,  we  also  show  the  results  from 
a  simple  block-diagonal  preconditioner  in  the  space  domain.  It  is  constructed  by 
inverting  the  block-diagonal  form  of  the  original  space-domain  moment  matrix.  The 
block  sizes  are  uniform  and  the  total  number  of  non-zero  elements  is  kept  the  same  as 
that  used  in  the  PWP  preconditioner.  We  find  that  the  iteration  number  grows  very 
rapidly  with  the  problem  size  if  no  preconditioning  is  applied.  With  the  PWP 
preconditioners  applied,  the  iteration  numbers  are  significantly  reduced  and  they  increase 
much  more  slowly  with  problem  size.  Furthermore,  the  PWP  preconditioner  performs 
better  than  the  equivalent  space  domain  preconditioner.  This  is  because  the  PWP 
preconditioner  captures  more  of  the  energy  in  the  original  moment  matrix.  Note  that  the 
CFBE  results  shown  in  Fig.  5(b)  have  smaller  iteration  numbers  than  the  results  from  the 
EFIE  cases  shown  in  Figs.  5(a)  and  5(c).  However,  the  effectiveness  of  the  PWP 
preconditioner  is  still  apparent.  In  Fig.  5(c),  we  also  plot  the  results  using  the 


conventional  wavelet  transform  (CWT)  basis  instead  of  the  PWP  basis.  The  CWT 
preconditioner  is  constructed  in  exactly  the  same  way  as  the  PWP  preconditioner,  except 
the  approximate  matrix  is  formed  from  the  elements  of  the  CWT  transformed  moment 
matrix  rather  than  the  PWP  transformed  one.  We  observe  that  the  PWP  preconditioner 
outperforms  the  CWT  preconditioner  since  the  PWP  basis  can  more  efficiently  represent 
the  moment  matrix.  Although  not  shown,  similar  comparisons  between  the  PWP  and 
CWT  preconditioners  are  also  found  for  Cases  1  and  2. 

Figs.  6(a)-(c)  show  the  convergence  behaviors  for  a  fixed  problem  size  (N=512) 
for  Cases  1-3,  respectively.  The  BiCGSTAB  algorithm  is  used  as  the  iterative  solver.  In 
each  figure,  the  residual  error  is  plotted  as  a  function  of  the  number  of  iterations  for  the 
moment  equation  without  any  preconditioning,  with  the  space  domain  and/or  the  CWT 
preconditioner,  and  with  the  PWP  preconditioner.  We  observe  that,  with  the  PWP 
preconditioning,  the  relative  residue  decreases  the  fastest.  In  addition,  the  convergence 
behaviors  are  more  stable.  Similar  convergence  behaviors  are  found  for  larger  problem 
sizes. 

Next,  we  investigate  the  effect  of  the  depth  index  K  in  the  wavelet  tree  on  the 
resulting  PWP  preconditioner.  For  a  scattering  problem  with  a  problem  size  of  N,  the 
maximum  depth  index  for  the  PWP  decomposition  tree  is  (logaN).  However,  this  does 
not  necessarily  lead  to  the  best  preconditioning  performance.  Fig.  7  plots  the  required 
iteration  number  versus  the  depth  index  K  used  in  constructing  the  PWP  preconditioner 
for  N=1024.  The  results  show  that  the  curves  flatten  considerably  after  K=5,  implying 
that  there  is  not  much  to  be  gained  beyond  this  depth  in  the  wavelet  tree.  The  reason  for 
this  phenomenon  is  because  although  the  maximum-depth  PWP  tree  generates  the 


sparsest  transformed  moment  matrix  [9],  the  long  PWP  basis  functions  resulting  from 
that  tree  shift  some  of  the  energy  into  the  non-diagonal  blocks,  making  the 
preconditioning  results  worse.  For  the  problem  sizes  we  have  studied  from  128  to  4096, 
the  optimal  PWP  tree  depths  are  usually  between  3  and  6.  Since  the  computational  cost 
of  constructing  the  preconditioner  is  proportional  to  the  maximum  depth  of  the  PWP 
decomposition  tree,  we  should  not  go  beyond  this  depth  when  constructing  the  PWP 
preconditioner. 

Finally,  we  examine  the  computational  cost  of  the  PWP  preconditioner.  Fig.  8 
shows  the  CPU  running  time  required  to  generate  the  PWP  preconditioning  matrix  [  Z  bd] 
for  Case  1  versus  problem  size  N.  The  direct  element  computation  in  equation  (10)  is 
carried  out  and  the  PWP  tree  depth  K  is  increased  slightly  from  3  to  6  as  the  problem  size 
is  changed  from  128  to  4096.  Fig.  9  shows  the  CPU  running  time  needed  to  implement 
the  multiple  matrix-vector  products  required  for  the  carrying  out  the  PWP 
preconditioning  in  (9).  We  exclude  the  time  needed  for  the  first  moment  matrix  and 
vector  product  in  (9),  as  that  operation  can  be  done  via  the  MLFMM.  We  find  that  the 
CPU  running  time  grows  at  rate  close  to  (9(NlogN)  in  the  both  cases.  This  is  consistent 
with  the  estimates  given  in  Section  IH  and  meets  the  objectives  we  have  set  forth  for  the 
preconditioner. 


V.  Conclusions 

A  preconditioner  for  the  moment  equation  based  on  the  pre-defined  wavelet 
packet  basis  has  been  proposed  in  this  paper.  Due  to  the  vanishing  moments  of  wavelet 
basis  functions,  the  PWP-based  moment  matrix  is  sparse  and  diagonally  concentrated. 


Consequently,  an  approximate-inverse  preconditioner  can  be  more  easily  designed  and 
constmcted  than  that  based  on  the  original  space-domain  moment  matrix.  The  PWP 
preconditioner  is  constructed  by  inverting  a  block-diagonal  form  of  the  PWP-based 
moment  matrix  and  transforming  the  resulting  matrix  back  into  the  space  domain.  It  has 
been  shown  that  the  complexity  for  the  construction  and  preconditioning  operation  can  be 
kept  under  O(NlogN)  in  both  computational  cost  and  memory  requirement.  Our 
numerical  results  showed  that  the  iteration  numbers  for  the  PWP-preconditioned  moment 
equations  are  significantly  smaller  and  grow  at  a  lower  rate  than  those  without 
preconditioning  or  preconditioned  using  a  space-domain  preconditioner.  In  addition,  the 
PWP  preconditioner  also  performed  better  than  the  equivalent  preconditioner  derived 
from  the  conventional  wavelet  basis.  We  should  point  out  that  in  this  paper  the  PWP 
algorithm  is  implemented  only  for  two-dimensional  problems.  Its  further  extension  to 
three-dimensional  (3-D)  problems  is  still  under  investigation,  with  the  goal  of  combining 
the  preconditioner  with  the  3-D  MLFMM  to  obtain  solutions  of  large-scale  problems. 

Finally,  although  we  have  chosen  the  block-diagonal  matrix  as  the  underlying 
structure  of  our  approximate-inverse  preconditioner,  we  have  found  that  a  block-banded 
matrix  [5]  or  a  threshold  transformed  matrix  [11,  20]  gave  better  preconditioning  results 
than  the  block-diagonal  form.  However,  the  inverses  of  those  matrices  are  dense, 
resulting  in  O(N^)  computational  bottlenecks  for  implementing  the  matrix-vector  product 
in  the  preconditioning  step.  A  worthwhile  topic  is  to  find  a  more  effective  approximate 
matrix  in  the  PWP  basis  domain  than  the  block-diagonal  one,  while  preserving  the 
D(NlogN)  total  computational  cost  and  memory  requirement. 
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Figure  1 .  (a)  Conventional  wavelet  decomposition  tree; 
(b)  Pre-defined  Wavelet  Packet  (PWP)  decomposition  tree. 


Figure  2.  The  pattern  used  to  construct  [  Z  bd]  from  [  Z  ]. 


(C) 


Figure  5.  Iteration  numbers  vs.  problem  sizes  for  solving  moment 
equations  with  different  preconditioning  methods  for  the  scattering 
structures  in  (a)  Case  1,  (b)  Case  2,  and  (c)  Case  3.  The  problem  sizes 
are  increased  by  proportionally  increasing  the  scatterer  sizes. 
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Figure  6.  Convergence  behaviors  of  the  preconditioned  systems 
versus  iteration  numbers  for  PWP  and  other  preconditioning  methods  in 
(a)  Case  1,  (b)  Case  2,  and  (c)  Case  3  with  N=512. 
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Figure  7.  Effect  of  the  PWP  decomposition  level  on 
iteration  number  for  Cases  1-3  using  the  PWP  preconditioner 

with  N=1024. 
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Figure  8.  CPU  running  time  vs.  problem  size  for  forming  the  PWP 
matrix  [  Z  bdl  using  ( 1 0)  for  Case  1 . 
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Figure  9.  CPU  running  time  vs.  problem  size  for  implementing  multiple 
matrix-vector  products  in  (9)  for  the  PWP  preconditioning  in  Case  1. 
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6  firom  a  single  mobile  user.  Equation  (1)  rquesents 
the  basic  form  of  an  ideal  steering  vector. 


au(0)= 


exp(-jk(x,  cos  6 + y ,  sin  9)) 
exp(“-jk(x2  COB  0 + y  2  sin  0)) 
cjq)(~jk(x  jCOsO+yjSin©)) 


cxp(-jk(x  cos  0 + y  M  sin  0)) 


(1) 


In  the  above  equation,  k  is  the  wavenumber  of  the 
incident  clectromagnedc  radiation  that  can  be 
expressed  as  2tiA..  Note  the  implicit  assumptions  of 
this  model  do  not  take  into  account  any  of  the  non¬ 
ideal  array  effects  given  in  Section  I.  Specifically,  the 
equal  gain  of  each  of  the  chamiels  shofws  diat 
amplitude  and  phase  mismatch  between  channels  is  not 
'considaed.  Also,  the  model  of  the  steeling  vector  in 
Equation  (1)  docs  not  include  any  terms  indicating  the 
retransmission  of  signal  oompiNients  fiom  one  array 
element  to  another  (i.e.  mutual  coupling).  Finally,  the 
model  assumes  that  the  clement  locations  (Xj,  yj  axe 
knofwn  exactly.  Iftbe  actual  antenna  element  locations 
are  pertnibed  by  a  significant  electrical  leagth,  the 
steering  vector  can  be  adversely  affected. 

A  distortion  matrix,  C(0)»  is  generally  used  to 
encapsulate  all  d  these  um-ideal  effects.  In  the 
literature,  this  matrix  is  estimated  from  collected 
measurement  data  [3,6]  and  is  applied  to  the  equation 
of  the  ideal  steering  vector  form  of  (1)  to  develop  a 
“compensatetT  steermg  vet^. 

ie(0)=C(0)5„(0)  (2) 


In  Section  HI,  we  describe  the  procedures  used  to 
solve  for  g  and  $0(6)  and  thus  detenninc  an 
estimate  for  die  actual  array  steeriiig  vector. 

HL  CAUBKATION  PROCESS 

A,  Cable  Measurements 

Figure  1  shows  a  diagram  illustrafing  how  th£ 
vector  g  is  determined.  There  are  two  sets  of  relative 
auqilitude  and  phase  infimnation  that  make  up  this 
detonnination.  The  first  set  corresponds  to  the  long 
mhling  leading  to  the  antenna  array  itsctf  on  the 
cdlular  tower  (rdfened  to  as  “array**  cables  in  Figure 
1).  The  secimd  ^  coriespoiids  to  the  combined  efGsct 
of  the  cables  used  to  connect  die  basestation  to  the 
array  c^Kes  (lefianred  to  as  “basestation”  cables  in 
Figure  1)  and  the  varying  impedance  seen  lookiiig  into 
the  ports  of  the  basestatiem  hardware. 


The  (Ustortion  matrix  is  generally  considerBd  to  be 
indqiendent  of  angle.  Thus,  mediods  iwing  this 
compensation  tedmique  do  libt  correct  for  tower 
effects  and  general^  assume  the  an^  dements  to  be 
very  similar,  if  not  idendcal. 

Mdhods  using  a  distortion  matrix  gcaerally 
consider  separately  the  issue  of  ami^itnde  and  phase 
mismatch  between  array  element  cabling.  To  take  this 
into  account,  we  ap|^  the  M-vector  g  that  contEdns 
the  complex  values  required  to  ccunpcasatc  tor  this 
mismatch  and  solve  tor  the  steering  vector 

(0)  (where  the  sttoscripi^i”  indicates  the  i***  vector 
oonqioDmit)  wdiich  can  be  used  in  DOA  analysis: 

aR.i(0)=giac.i(6)  (3) 


It  is  leladvely  BtraighCtorwazd  to  determine  the 
relsfive  ftmpiifndft  and  idiase  response  of  the  array 
cables  using  a  netwodc  analyzer.  We  refer  to  this 

VBCtoras  T  and  it  only  needs  to  be  detomined  once  fm 
cadi  set  of  array  otolesL  This  vector  conlains  the 
relative  amplitude  and  phase  reqxmse  of  each  of  the 
cables  with  reflect  to  a  given  reference  cable  (that  has 

a  1Z0°  as  its  correqxmding  entry  in  the  I  vector). 

It  is  slightly  mme  difficult  to  determine  the 
idafivo  amplitnde  and  phase  response  of  the 
basestation  cables  and  RF  impedance  lookiiig  into  the 
basestation  Q’^rred  to  as  the  S  vector)  since  this 
measnrement  must  be  made  every  time  the  basestation 
is  powered  tip.  Again,  rh«  nmasuienieat  is  made  with 
lesped  to  a  reference  dumncL  Using  a  power  splitfer 
with  a  known  relative  amplitude  and  phase  nspoaso, 
p ,  (deiBnnined  via  network  analyzer)  an  RF  signal 
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Crom  a  signal  generator  is  injected  into  the  basestation 
cables.  The  resulting  data  returned  by  the  basestadon 
normalized  by  the  reference  channel  value  are  the 
elements  of  the  i  vector,  g  is  then  determined  by: 


Pi 

The  right  hand  side  of  Equation  (4)  contains  the 
relative  amplitude  and  phase  mismatch  of  the  array 
cables  and  basestation  cibles/haixtware.  The  e£fect  of 
the  power  splitter  most  be  removed  since  the  splitter 
was  included  in  the  measurement  of  the  s  vector  but  is 
absent  from  normal  system  operation. 


B.  CEM  Simuiations 


There  are  two  different  ways  Id  solve  for  * 
The  first  is  to  estimate  the  C  matrix  algorithmically  in 
Equation  (2)  as  is  done  in  [3]  and  then  use  this  estimale 
to  determine  i^(9).  An  aherriative  inethod  [7]  iises 
CEM  simtilationa  [8]  to  determine  ic(B)  directly. 


The  CEM  padcage  used  in  this  study  is  NEC 
(Numerical  Electromagnetics  Code),  which  was- 
originally  developed  at  the  Lawrence  Livermore 
Laboratory  [9]  to  perform  moment  method  analysis  to 
modd  the  interaction  between  electromagiietic  fields 

wire  segments.  A  wise  segment  model  of  an 
arbitrary  geometry  antenna  array  can  be  specified  in  an 
input  file  to  NEC.  Rgnrc  2  shows  the  model  of  a 
dnglft  antenna  array  dement  used  in  simulations.  An 
arr^  of  these  elements  are  placed  according  to  array 
geometry,  incident  signals  arc  launched,  one  at  a 
time,  firom  ail  dircctkms.  The  air^  output  vector  for 
each  of  these  incident  signals  is  (0) . 

IV.  E3CPERIMENT  SETUP 

The  array  considered  in  this  study  is  a  uniform 
Cfrc"!**'  antenna  array  (UCA)  operating  at  1.8  GHz. 
Field  measurements  were  made  at  the  Pkkle  Researdi 
Campus  at  the  University  of  Texa^  at  Austin.  To 
reduce  multipath  signal  Gomponeiits,  measurements 
were  njadc  in  a  field  with  buildings  for  off  in  the 
distance.  An  RF  signal  generator  attached  to  a  dipole 
antenna  was  used  to  represent  a  mobile  user. 
Measunmients  were  made  with  this  mobile  user  at 
several  different  locations  relative  to  the  basestation. 
DOA  analysis  of  the  received  data  is  made  with  the 
Multiple  Signal  Classification  ^fUSIQ  algorithm 
[10]. 

V.  RESULTS 

A,  Pre-CaJibration  J^patial  Sji>ectntm 


Figure  3  -  MUSIC  Spatial  spectrum -No  Calibration 

Figure  3  shows  the  MUSIC  spsiM  spectrum  due 
to  a  mobile  user  located  at  idative  to  the 
baaestatkm.  As  seen  fay  the  above  plot,  tiie  lobes  in  the 
iKtfinalized  qiatiai  qKctrum  do  rid  omte^wnd  at  all  to 
the  signal  of  the  desired  mobile  user.  In  foct,  the  lobes 
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in  this  spatial  spectnmi  do  not  correspond  to  any 
known  signal  or  multipath  component,  which  is  not 
surprising  due  to  the  lade  of  any  kind  array 
calibration. 

B.  Amplitude  and Phase^tsmatch  Compensalion 

Figure  4  shows  the  MUSIC  spatial  spectrum  of 
Figure  3,  with  compensation  added  for  array  cables 
and  basestaiion  cd^les/hardware,  as  discussed  in 
Section  nL  A.  As  this  figure  shows,  the  mobile  user  at 
80°  can  now  be  restdved  easily.  However,  again  there 
is  a  lobe  at  280°  that  does  not  correspond  to  ar^ 
known  signal  component  or  multipath. 


Figure  4  -  MUSIC  Spatial  Spectrum  -  Cable 
Cah1)ration 


Figure  5  -MUSIC  Spatial  plectrum  -Full  Am^ 
Calibration  Applied 

C  Mutual  Couplu^  Compensation 

Figure  5  displi^  the  resnh  of  adding  the  mutual 
oocqding  compensation  discussed  in  Section  in.B.  to 
the  ^lectnim  shown  in  Figure  4.  Specifically,  the 
results  in  this  figure  includes  calibration  for  amiditiide 


and  phase  mismatch  due  to  array  cd)le8  and 
basestation  cables/hardware,  as  well  as  calibration 
including  mutual  coupling  effects.  While  there  is  still 
an  extraneous  ld}e  in  this  Figure,  it  is  not  as  large  as 
the  lobes  in  Figures  3  or  4,  which  could  easily  be 
counted  as  multipath  signal  components  fitom  diffisient 
DOAs.  In  addition,  comparison  of  Figure  4  and  5 
shows  a  narrower  peak  in  Figure  5  that  corresponds  to 
the  DOA  of  the  desired  mobile  user.  This  general 
result  was  observed  in  all  tested  cases. 

VL  CONCLUSIONS 

The  results  of  this  paper  quantify  the  effects  of  a 
multi-step  sm^  antmma  calibration  process.  As 
shown  by  noting  the  relative  improvement  from 
Figures  3  and  4  and  the  improvement  from  Fignies  4 
and  5,  it  is  dear  that  anqilitude  and  phase  mismatch 
effects  in  array  cables  and  basestation  cables/hardware 
are  much  more  significant  than  mutual  coupling 
effects.  This  is  not  to  imply  that  mutual  emoting 
calibration  is  not  necessary  since  Figure  5  cleariy 
shows  that  lobes  that  could  have  been  mistaken  for 
multipath  signal  energy  were  effectively  reduced.  In 
addition,  mutual  coupling  cafibration  makes  the  peak 
in  the  spatial  spectrum  corresponding  to  the  desired 
mobile  user  narrower. 

The  procedure  presented  in  this  p^per  is  important 
because  it  does  not  require  expensive  and  involved 
site-^Kcific  calibiatifHi.  The  network  analyzer 
measurements  can  all  be  made  relatively  easily  and  the 
CEM  simuiations  can  be  performed  ofr-Une  using  the 
free  and  well-accqitBd  package  NEC.  The 
combination  of  these  two  steps  allow  for  very  effective 
arr^  calibration. 

ACKNOWLEDGEMENTS 

This  work  is  siqipoited  by  the  Texas  Hi^ier 
Education  Coordinating  Board  under  the  Texas 
Advanced  TechnoU^  Program  and  by  the  Office  of 
Naval  Research  under  contract  no.  N00014-98-1-0178. 

REFERENCES 

[1]  Y.  Chen,  A  Chsmg,  and  R  Lee,  ‘‘Array 
Calibration  Methods  for  Sensor  Position  and 
Pointing  .  Erron,**  Microwave  and  Optical 
Technology  Letters,  \o\.  26,  pp.  132-137, 2000. 

[2]  B.  Kang,  R  Subbaram,  and  B.  Steinberg, 
‘Tmpioved  Adaptivo-Beamfoimin^  Target  fm 
Self-Calibrating  aDistorted  Phased  Array,**  IEEE 
Transactions  on  Antennas  and  Propagation,  vol. 
38,  pp,  186-194, 1990. 


PI  C.  See,  “Sensor  array  calibration  in  the  presence 
of  mutual  coupling  and  unknown  sensor  gains 
and  phased”  Electronics  Letters,  voL  30,.  pp. 
373-374,  1994. 

[4]  R.  Ertttl,  Z.  Hu,  and  J.  Reed,  “Antenna  Array 
Hardware  Amplitude  and  Phase  Compensation 
Using  Basd>and  Antenna  Array  Outpi^”  1999 
IEEE  Vehicular  technology  Conference 
Proceedings,  voL  3,  pp.  1759-1763, 1999. 

|5]  G.  Tsoulos  and  M.  Beadi,  Xalibradon  and 
Lineality  issues  ibr  an  Adaptive  Antenna 
System,”  1997  IEEE  Vehicular  Technology 
(inference  Proceedings,  vd.  3,  pp.  1597-1600, 
1997. 

[6]  A  Lemma,  E.  Deprettere,  and  A.  van  der  Veen, 
“Experimental  Analysis  of  Antenna  Coupling  for 
High-Resolution  DOA  Estiination  Algorithms,” 
1999  IEEE  T*  Workshop  on  Signal  Processing 
Advances  in  Wireless  Communications,  pp.  362- 
365,  1999. 

[7]  T.  Su,  K.  Dandekar,  and  R  ling,  “Simulation  of 

Mutual  Coupling  Effect  in  Circular  Arrays  for 
Direction  Finding  Applications,”  hPicrawtxvt  and 
Optical  Technology  Letters,  vol.  26,  331-336, 

2000, 

[5]  K.  Pasala  and  E.  Friel,  “Mutual  Coupling  Ef&cts 
and  Their  Reduction  in  Wkkband  Directioa  of 
Arrivai  Estimation,”  IEEE  Thmsactions  on 
Aerospace  and  Electronic  Systems,  vol.  30,  pp. 
1116-1122. 

[9]  NEC-2  Manual,  Lawrence  liveniiore  National 
Laboiatoty,  1996. 

(101  R  Krim  and  M.  \^bcrg,  “Two  Decades  of  Array 
Signal  Processing  Research:  The  Parametric 
Approach,”  IEEE  Signal  Proces^ng  kfagazine, 
pp.  67-94,  1996. 


Paper  [26] 


Design  of  Dual-Band  Microstrip  Antennas  Using  the  Genetic  Algorithm 


H.  Choo,  and  H.  Ling 

Department  of  Electrical  and  Computer  Engineering 
The  University  of  Texas  at  Austin 
Austin,  TX  78712-1024  U.S.A 
E-mail:  hschoo@ece.utexas.edu 


Abstract: 

We  report  on  the  use  of  the  genetic  algorithm  (GA)  to  design  patch  shapes  of  microstrip  antennas  for 
dual-band  applications.  We  employ  a  full-wave  electromagnetic  solver  to  predict  the  performance  of 
microstrip  antennas  with  arbitrary  patch  shapes.  Two-dimensional  chromosome  is  used  to  encode  each 
patch  shape  into  a  binary  map.  GA  with  two-point  crossover  and  geometrical  filtering  is  implemented 
to  achieve  efficient  optimization.  The  GA-optimized  designs  are  built  on  FR-4  substrate  and 
measurement  results  show  good  agreement  with  the  numerical  prediction.  The  optimized  patch  design 
achieves  good  impedance  match  at  both  frequencies.  The  patch  shape  is  further  optimized  to  broaden 
the  bandwidth  at  one  of  the  frequencies.  It  is  also  shown  that  a  frequency  ratio  between  the  two 
frequency  bands  ranging  from  1.2  to  2  can  be  achieved  using  the  GA  design  method. 


Introduction: 

With  the  growing  demand  of  wireless  applications,  microstrip  antennas  that  can  operate  in  two  or  more 
frequency  bands  are  in  demand.  Various  dual-band  methods  for  microstrip  antennas  have  been 
proposed  to  date.  For  example,  multi-layered  structures,  parasitic  patches  and  shorting  posts  are  some 
of  the  well-known  techniques  for  achieving  dual-band  operation  [1].  In  this  paper,  we  examine  the  use 
of  genetic  algorithms  (GA)  to  design  optimal  shapes  for  microstrip  antennas  to  achieve  dual-band 
operation.  The  attractiveness  of  GA  shape  optimization  is  that  dual-band  performance  can  be  achieved 
with  little  increase  in  overall  volume  or  manufacturing  cost. 

The  design  of  dual-band  microstrip  antennas  using  GA  was  first  addressed  by  Johnson  and  Rahmat- 
Samii  [2].  Air  substrate  was  used  in  their  study.  In  this  work,  we  focus  on  FR-4  as  the  substrate 
material,  since  it  is  the  most  commonly  used  material  in  wireless  devices.  In  our  GA  implementation,  a 
two-point  crossover  scheme  involving  three  chromosomes  is  used.  Furthermore,  geometrical  filtering 
is  applied  to  each  chromosome  to  obtain  a  more  realizable  shape.  The  GA-optimized  shapes  for  dual¬ 
band  operation  are  presented.  It  is  also  shown  that  arbitrary  frequency  ratios  between  the  two 
frequency  bands  ranging  from  1.2  to  2  can  be  achieved  through  the  GA  design. 


GA  Optimization: 


GA  is  implemented  to  optimize  the  microstrip  patch  shape  in  order  to  achieve  dual-band  operation.  In 
our  GA,  we  use  a  two-dimensional  (2-D)  chromosome  to  encode  each  patch  shape  into  a  binary  map 
[3].  The  metallic  sub-patches  are  represented  by  ones  and  the  no-metal  areas  are  represented  by  zeros. 
Since  it  is  more  desirable  to  obtain  optimized  patch  shapes  that  are  well  connected  from  the 
manufacturing  point  of  view,  a  2-D  median  filter  is  applied  to  the  chromosomes  to  create  a  more 
realizable  population  at  each  generation  of  the  GA. 

To  evaluate  the  performance  of  each  patch  shape,  we  use  a  fiill-wave  periodic  patch  code  adapted  from 
a  frequency  selective  surface  code  [4].  The  formulation  of  the  code  is  similar  to  that  described  in  [5]. 
The  electromagnetic  analysis  is  carried  out  by  using  the  electric-field  integral  equation  (EFIE).  The 
periodic  Green’s  function  for  layered  medium  is  the  kernel  of  the  integral  equation.  Rooftop  basis 
functions  are  used  to  expand  the  unknown  current  on  the  metal  patch  and  fast  Fourier  transform  (FFT) 
is  used  to  accelerate  the  computation  of  the  matrix  elements.  To  reduce  the  matrix  fill-time,  matrix 
element  calculation  is  done  only  once  and  stored  before  the  GA  process.  Because  of  the  assumed 
periodicity  in  this  patch  code,  we  use  a  period  that  is  greater  than  one  wavelength  to  avoid  coupling 
between  the  adjacent  patches  for  single  patch  simulation. 

The  design  goal  is  to  produce  good  antenna  radiation  at  two  frequency  bands  by  changing  the  patch 
shape.  The  low  frequency  is  chosen  at  1.9  GHz  and  the  higher  one  is  chosen  at  2.85  GHz.  To  achieve 
dual-band  design,  an  additive  cost  function  is  defined  as 
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+  if  St^(dB)>- 10 dB 

0  if  Sn(dB)< -10 dB 
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The  first  part  of  the  cost  function  accounts  for  the  impedance  mismatch  and  is  defined  as  the  average 
of  those  return  loss  (Sn)  values  that  exceeds  -lOdB  (i.e.,  VSWR=2:1)  within  the  two  frequency  bands 
of  interest.  The  second  part  of  the  cost  function  accounts  for  the  total  metal  loss  (dB)  generated  by  the 
current  flowing  on  the  patch.  The  conductivity  of  aluminum  (a  =  3.82e7  S/m)  is  used,  as  the 
microstrip  in  our  measurements  were  built  using  aluminum  tape. 


Based  on  the  cost  function,  the  next  generation  is  created  by  a  reproduction  process  that  involves 
crossover,  mutation  and  2-D  median  filtering.  A  two-point  crossover  scheme  involving  three 
chromosomes  is  used.  The  process  selects  three  chromosomes  as  parents  and  divides  each 
chromosome  into  three  parts.  Intermingling  the  three  parent  chromosomes  then  makes  three  child 
chromosomes.  It  is  found  that  the  two-point  scheme  exhibits  a  faster  convergence  behavior  than  the 


one-point  scheme.  The  GA  process  is  iterated  until  the  cost  function  converges  to  a  minimum  value. 

Results: 

Fig.  1(a)  shows  the  GA-optimized  microstrip  shape  for  dual-band  operation.  A  72mm  x  72mm  square 
design  area  in  which  the  metallic  patch  can  reside  is  discretized  into  a  32  x  32  grid  for  the 
chromosome  definition.  The  thickness  of  the  FR-4  substrate  (dielectric  constant  of  4.3)  is  1 .6  mm.  In 
the  GA-optimized  shape,  the  dark  pixels  are  metal  and  the  white  pixels  have  no  metal.  The  white  dot 
shows  the  position  of  the  probe  feed.  Fig.  1(b)  shows  the  predicted  return  loss  (Sn  in  dB)  of  the 
resulting  microstrip  antenna.  Good  matches  are  exhibited  by  the  return  loss  curve  at  the  design 
frequencies  of  1.9  GHz  and  2.85  GHz.  The  bandwidths  at  the  two  design  frequencies  are  about  4% 
and  1.4%,  respectively. 

Previously,  we  have  used  GA  to  achieve  broadband  operation  for  a  single-band  microstrip  [6].  The 
resulting  bandwidth  on  1.6  mm  FR-4  substrate  is  about  8%.  Here,  we  try  to  increase  the  bandwidth  of 
the  dual-band  microstrip  antenna  from  the  above  design.  The  low  frequency  (1.9  GHz)  is  chosen  to  be 
our  target  for  broadbanding,  while  the  bandwidth  of  high  frequency  is  kept  the  same.  During  the  GA 
iterations,  we  gradually  increase  the  frequency  range  centered  at  1.9  GHz  in  our  cost  function 
definition  until  the  broadband  design  is  achieved.  Fig.  2(a)  shows  a  bandwidth-enhanced  dual-band 
result.  To  experimentally  verify  the  GA  design,  we  have  built  and  measured  such  a  microstrip  patch. 
Fig.  2(b)  shows  the  measurement  and  simulation  return  loss  for  the  shape  in  Fig.  2(a).  Good  agreement 
is  observed  between  the  measurement  and  simulation  results  (the  square  dots  show  the  predicted 
values).  Figs.  2(c)  and  (d)  show  the  measured  boresight  radiations  (S21  dB)  for  the  two  diagonally 
oriented  polarizations.  Also  plotted  in  the  insets  are  the  predicted  current  distributions  at  three 
frequencies  of  interest.  We  note  that  near  1.9  GHz,  there  are  two  modes  with  orthogonal  current 
directions  at  two  closely  spaced  frequencies,  leading  to  the  broadening  of  the  impedance  bandwidth  in 
Fig.  2(b).  At  2.9  GHz,  only  a  single  mode  exists. 

In  the  above  examples,  the  ratio  between  the  two  frequency  bands  is  chosen  to  be  1:1.5.  We  next  carry 
out  the  dual-band  design  steps  using  GA  to  achieve  different  frequency  ratios  between  the  low  and  the 
high  frequency  band.  During  the  design,  we  keep  the  low  frequency  band  fixed  around  1.9  GHz, 
while  varying  the  design  frequency  for  the  high  band.  Frequency  ratios  from  1:1.2  to  1:2  are 
attempted  using  GA.  Figs.  3(a),  (b)  and  (c)  show  the  optimized  shapes  and  the  corresponding  return 
loss  versus  frequency  curves  for  the  frequency  ratios  1:1.2,  1:1.4  and  1:1.74,  respectively.  We  find 
that  it  is  possible  to  use  the  GA  approach  to  cover  the  entire  dual-band  frequency  ratio  from  1 : 1 .2  to 
1:2. 

Conclusions: 

Optimized  patch  shapes  for  dual-band  microstrip  antennas  on  thin  FR-4  substrate  have  been 
investigated  using  the  genetic  algorithm.  The  resulting  microstrip  designs  show  good  operation 


characteristics  for  dual  frequencies.  The  bandwidth  of  the  microstrip  can  be  further  broadened  by 
optimization.  This  result  has  been  verified  by  laboratory  measurement.  Finally,  it  has  been  shown  that 
the  frequency  ratio  between  the  two  bands  ranging  from  1 : 1 .2  to  1:2  can  be  achieved  using  the  same 
GA  methodology. 
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Fig.  1.  GA  optimized  dual- band  microstrip  antenna  (frequency  ratio  1:1.5). 
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IntrodBClioii 

Genetic  algorithm  (GA)  has  been  applied  to  the  design  of  imilri-layeiEd  planar 
and  cylindrical  absorber  stntcturcs  [11.  Comigaicd  coatings  with  oon^planar  shape 
protile  offer  an  additicmal  degree  of  design  freedom  and  have  been  analyzed  in  PI.  fa 
ibis  papier,  we  apply  GA  nr  design  (be  optimal  shape  for  a  corrugated  coating  under 
obtique  incidence  (see  Fig.  1).  To  predict  the  perfoniunoe  of  each  shape,  wc  employ  a 
fol}*wave  electromagnetic  solver.  GA  with  two^point  crossover  and  gcometiical  filtering 
is  implememed  to  achieve  effective  optimization.  The  optimized  designs  for  different 
polarizaiions  are  ftfcsented. 

G  A  Optimization 

GA  is  impleniented  to  optimize  the  shape  of  a  corrugated  absorber  under  oblique 
incidenceL  The  absorber  consists  of  a  sin^  layer  of  lossy  material,  h  has  one- 
dimensional  periodicity  along  the  X'4imeosion  and  is  invariant  along  the  ZKliincnsion 
(Fig.  I).  The  incidence  wave  direetion  is  deftaed  by  Gd  and  The  elevation  angle 
6ld  is  measured  from  grazing  while  ^  azimuth  ai^le  measured  from  the  zraxis. 

To  encode  each  possible  ahsotber  shape  into  a  chrornosonie,  the  height  of  the  coating  at 
each  point  along  die  x-direction  U  lepnacnt  as  a  bmaiy  string.  To  achieve  more 
realizable  shape,  wc  apply  a  sliding  wii^w  filter  that  operates  like  a  bw-poss  spatial 
tiller  during  each  generation  of  the  GA.  fa  addition,  we  enforce  f^nninttry  con^raint  on 
the  shape. 

To  evaluate  the  performance  of  each  absorber  shape,  wc  use  a  full-wave 
electromagnetic  simulation  code  based  on  the  boundary-integral  formulation  [2].  The 
design  goal  is  to  choose  the  coating  protile  that  gives  rise  to  the  best  absorbing 
characteristics  over  the  frequency  band  of  interest  The  design  frequency  band  is  chosen 
to  be  from  8  GHz  to  18  GHz,  and  the  maximum  height  of  the  grating  is  restricted  to 
8mm.  Associated  with  this  design  goal,  we  deflne  the  cost  function  as: 
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^  ^  „  fr(^B)+2(WB 

•  Cost  —  — 2  where  Pn-< 

(0 

Based  on  the  cost  function,  the  next  generation  is  created  by  a  Tcproduction  process  that 
involves  crossover,  mutation  and  sliding  window  fihering.  A  two-point  crossover 
scheme  involving  three  chromosomes  is  used-  The  process  selects  three  chromosomes  as 
parents  and  divides  each  chromosome  into  three  parts.  Intermingling  the  three  parent 
chromosomes  then  makes  three  child  chromosomes.  The  reproduction  process  is  iterated 
until  the  cost  function  converges  to  a  minimum  value. 

Results 

Standard  magram  is  used  for  the  coating  mateiial-  We  consider  the  case  when  the 
incident  angles  are  set  to  and  The  bottom  of  the  coating  is  backed  by  a 

conducting  ground  plane.  The  period  of  the  profile  is  set  to  2.032mm  and  is  discreti^ 
into  32  points.  The  height  of  the  groove  at  each  point  is  described  by  a  6>bit  number  that 
ranges  between  0  and  8mm  (in  64  steps).  First,  we  consider  the  case  when  only  the  W 
reflection  coefficient  is  used  in  the  cost  function  definition.  Fig.  2(a)  shows  the  resulting 
GA-optimized  shape,  which  closely  resembles  the  triangular  sh^  Fig.  2(b)  is  plot  of 
the  simulated  reflection  coefTtcients  (in  dB)  versus  frequency  for  the  optimized  shape. 
We  see  that  the  reflection  coefficient  of  the  W  polarization  nearly  meets  the  -20  dB 
design  goal  over  the  entire  frequency  band  from  8GHz  to  1 8GHz.  The  HH  polarization  is 
not  optimized  and  shows  a  much  higher  reflection  coefficienL  Next,  we  consider  the 
reverse  situation  when  only  the  HH  reflection  coefficient  is  used  in  the  cost  function. 
Hg.  3(a)  shows  the  resulting  GA-optimized  shape.  It  is  noted  that  the  optimal  shape  of 
ovrugation  resembles  a  rectangular  profile.  Fig.  3(b)  shows  the  siimilated  reflection 
coeffictenu  (in  dB)  versus  frequency  for  the  optimized  sh^  for  both  the  W  and  the  HH 
polarization.  The  leflectioa  coefficient  of  the  HH  polarization  is  less  than  -15  dB  for 
nearly  the  entire  frequency  band  of  interest  In  the  third  example,  we  optimize  the  shape 
by  using  the  average  of  tte  reflection  coefficients  from  the  HH  and  VV  polarizations  in 
the  cost  function.  The  resulting  shape  is  shown  in  Fig.  4<a).  As  we  have  seen  from  the 
last  two  examples,  the  design  for  the  HH  polarization  is  harder  than  that  for  the  VV 
polarization.  Therefore  in  this  case,  the  cost  is  dominated  by  the  HH  consideration,  and 
the  resulting  GA^optinuzed  shape  is  not  that  different  from  the  optimized  shape  for  the 
HH-polarization  shown  in  Rg.  3(a).  As  a  final  example,  we  impose  an  additional 
constraint  such  that  the  complemejitary  air  region  has  the  same  shape  as  the  absorber. 
Prcsunnably,  this  means  wc  can  make  two  usable  pieces  of  absorber  by  cutting  them  from 
a  single  thick  piece.  The  optimized  shape  and  the  reflection  coefficient  are  presented  in 
Figs.  5(a)  and  5(b),  respectively. 


i/r(dB)2:-20d» 

i/r(dB)<-2(WB 


GA<optifnized  shapes  for  a  corrugated  absorber  under  oblique  incidence  have 
been  presented  for  diiTercot  incident  polarizations.  The  designed  absorber  shape  for  the 
W  p^arization  closely  resembles  a  triangular  profile,  while  that  for  the  HH  polarization 
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resembles  t  rectangular  profile.  An  optimized  absorber  profile  with  an  additional 
compleinentary  constraint  to  achieve  a  more  efficient  material  usage  is  also  reported. 
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fig.  2  (a)  GA  optimized  shape  for  the  VV  pol.  (b)  Reflection  coeffidcat  (dB)  versus  ftequcocy 


710 


(»)  . 


Period  =  2.032inm 


(b) 


Frequency  (Ghz) 


Fig.  3  (a)  GA  optimized  shape  for  the  HH  pol.  (b)  Reflection  coefficient  (dB)  versus  frequency 


(a) 


Li 


Period  «  2.032  nom 


(b) 


Fig.  4  (a)  GA  optimized  shape  for  the  average  of  VV  and  HH  pol.  (b)  Reflection  coefficient  (dB) 
versus  frequency 
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Fig.  5  (a)  GA  oprimized  shape  with  the  complementary  constraint,  (b)  Reflection  coefficient  (dB) 
versus  frequency 
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Abstract  -  The  application  of  Genetic  Algorithm  (GA)  optimization  to  the  design  of  ultra  ^wideband 
antennas  for  use  in  midtiservice  wireless  devices  is  de^ribed.  The  antennas  explored  in  this  study  are 
monopoles  with  planar,  Jully  -  metal  radiating  elements.  The  key  design  constraints  and  criteria  to  be 
considered  in  the  development  of such  antennas  are  discussed  Presented  and  discussed  are  results  from 
an  application  of  a  GA  optimizer  to  the  design  of bow  -  tie  and  reverse  bow-tie  monopole  antennas. 


1.  INTRODUCTION 

The  trend  in  the  wireless  industry  is  towards 
enabling  mobile  end-terminal  devices  to  achieve 
multi-service  performance.  To  date  this  has 
included  access  to  both  analog  and  digital  cellular 
bands  between  0.8  -  0.9  GHz  and  L8  -  1.9  GHz.  It 
will  be  necessary,  however,  to  operate 
simultaneously  at  other  bands  such  as  1.8  to  2.2 
GHz  for  3G  services  and  2.4  GHz  to  allow  access 
to  WLAN  and  technologies  such  as  Bluetooth. 
Positioning  information  is  also  desirable  requiring 
access  to  the  GPS  bands  at  1.227  GHz  and  L575 
GHz.  These  devices  in  the  future  may  also 
implement  some  form  of  the  emerging  Impulse 
Ultra-Wideband  (lUWB)  radio  for  supplemental, 
short  range  communications  and  high-rcsolution 
ranging  or  positioning  [I].  lUWB  ofiers  distinct 
advantages  for  reduced  power  consumption,  design 
complexity,  and  size.  Ibis  technology,  which  uses 
an  impulse-based  carrier,  requires  an  ultra-wide 
instantaneous  bandwidth  on  ^  order  of  60%  to 
100%. 

Clearly  it  desirable  to  design  the  smallest  possible 
end-terminal  to  access  all  of  these  services 
simultaneously.  The  antenna  imposes  a  limitation 
on  the  miniaturization  that  is  achievable  in  such  a 
device.  Ideally,  a  single,  i^ysically  small  antenna 
is  desired  that  can  operate  over  all  of  the  frequency 
bands  listed  above.  This  requires  that  the  antenna 
function  at  cadi  of  the  bands  between  0.8  and  2.4 
GHz  as  well  as  over  the  entire  band  associated  with 
lUWB.  Good  impedance  matdiing  and  high 
radiation  effideiicy  roust  be  achieved  over  tiiis 
entire  bandwidth  In  otder  to  guarantee  satisfoctory 
signal  receptioD  and  maintain  low  power  operation. 
lUWB  service  also  requires  a  linear  phase  response 
for  efiflcient  operation.  These  requhmencs  impose 
severe  restrictions  on  the  types  of  antennas  that 
may  used  for  such  a  multi-service  wireless  device, 


especially  when  small  size  is  of  key  concern. 

Theie  have  been  a  wide  variety  of  antennas 
proposed  for  multi-service  wireless  applications. 
For  instance,  a  Planar  Inverted  F  Anteima  (PEFA) 
patch-type  design  cqiablc  of  covering  the  0.9,  1.8, 
and  2.4  GHz  bands  has  been  proposed  [2].  While 
patch  designs  offer  small  and  /  or  conformal 
solutions,  they  are  all  inherently  band-limited. 
Monopoic  designs  using  frequency  independent 
elements  offer  a  number  of  wide  operating  bands 
[3].  Other  frequency  independent  designs  such  as 
spirals  offer  “flat"  instantaneous  matching 
bandwidths  as  wide  as  5:1  [4].  Frequency 
independent  antennas,  however,  are  inherently 
dispersive,  which  does  not  satisfy  the  phase 
linearity  requirement  of  lUWB.  Recent  studies 
indicate  that  planar,  fiilly-metal  monopoles 
(PFMM)  may  offer  very  vtide  instantaneous 
bandwidths  (as  much  as  12:1)  with  a  relatively 
small  size.  As  their  operating  principle  is  much 
like  that  of  the  traditional  wire  roonopole,  these 
elements  should  also  demonstrate  good  phase 
linearity. 

The  most  common  PFMM  antenna  is  the  bow-tie 
(BT)  which  has  long  been  known  to  possess 
broadband  performance.  The  BT  antenna  is  the 
planar  fonn  of  the  conical  anteima  for  which  a 
closed-form  analysis  is  possible.  Its  perfbnnance 
can  therefore  be  understood  in  terms  of  that  of  the 
cone.  Also  the  performance  of  Unite-sized  BT 
antennas,  both  in  terms  of  input  impedance  and 
radiation  patterns,  has  been  characterized 
experimeinally  m  the  literature  [5].  From  this  data, 
it  can  be  determined  tiiat  there  are  two  key 
parameters  in  the  design  of  BT  antennas,  the 
dement  height,  and  tiie  element  flare  angle.  The 
hei£^  essentially  deteimmes  the  operating  mode 
and  the  lower  frequency  limit  of  flie  antenna,  while 
the  flare  angle  control  the  variation  of  the  input 
impedance  over  frequenqr,  the  high  frequency 
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impedance  value,  as  well  as  the  resonance 
bandwidth.  From  this  analysis,  it  can  be  seen  that 
the  design  of  BT  antennas  is  rather  straightforward. 
For  this  reason,  BT  antennas  have  been  used 
successfully  in  many  practical  broad-band  radio 
systems. 

Other  PFMM  element  shapes  such  as  squares  and 
circles  have  proposed  for  broadband  operation  [6]. 
These  designs  differ  from  the  BT  in  that  they  must 
be  raised  above  the  ground  plane  by  a  fixed  height 
in  order  to  avoid  shorting  to  it.  Their  performance 
has  yet  to  be  fully  explored  due  to  die  lack  of  an 
analytical  model,  and  the  only  recent  emergence  of 
experimental  data.  For  these  reasons,  there  exists 
no  clear  design  methodology  for  such  element 
shapes,  and  no  means  by  which  to  c^thnize  their 
performance. 

Therefore,  a  study  was  proposed  to  evaluate  these 
new  shapes  using  the  Genetic  Algorithm  (GA) 
approach.  The  GA  was  chosen  due  to  its  ability  to 
optimize  problems  involving  multiple  design 
parameters  as  well  as  conn^lex,  and  conflicting 
design  criteria,  and  to  converge  to  a  nearly  global 
solution  in  doing  so.  There  has  also  been  much 
success  in  recent  years  in  applying  the  GA  to 
antenna  design  [7].  In  this  study,  the  GA  is  used  to 
determine  the  optimal  dimensions  of  the  selected 
element  shape  in  order  to  fulfill  the  given 
bandwidth  requirement  By  repeating  this 
procedure  for  increasing  bandwidth  requirements, 
it  was  hoped  that  a  design  methodology  for  the 
clement  would  emerge  and  die  absolute  limits  of  its 
broadband  performance  could  be  determined. 

n.  APPROACH 

The  element  chosen  for  initial  study  was  the 
reverse  bow-tie  (RBT).  Tliis  shape  was  chosen  due 
to  its  similarity  to  the  traditional  BT,  which  would 
allow  a  direct  comparison  of  the  performance  of 
identically  shed  elements.  Figure  1  compares  the 
two  configurations.  As  can  be  seen,  the  clement 
height,  H,  and  flare  angle,  oc,  are  design  parameters 
for  both  elements,  while  the  RBT  has  the  additional 
parameter  of  feed  height,  hf. 

A  GA  optimization  algorithm  was  written  that 
operates  on  encoded  combinations  of  these  design 
parameters.  Using  a  numerical  electromagnetic 
simulation  code,  the  response  of  individual  antenna 
designs  refuesented  by  these  parameter 
combinadons  is  determined.  From  this  data,  a  cost 
function  is  defined  which  assesses  the  relative 
goodness  of  each  design  according  to  some  pre¬ 


determined  criteria.  New  designs  are  generated  by 
mixing  parameter  information  between  different 
designs  through  die  crossover  operator,  and 
randomly  changing  parameter  information  through 
the  mutation  operator.  The  performance  is  again 
evaluated  for  these  new  designs.  This  process  is 
repeated  until  the  original  design  goal  is  met 

a 


a.  b. 

FIGURE  1.  Monopolc  antenna  dements  under  study. 

a.  bowtie 

b.  revemebowde 

This  same  algoridim  was  applied  to  botii  the  BT 
and  RBT.  To  maintain  consistency,  the  design 
parameters  used  for  both  designs  were  height,  flare 
angle,  and  feed  height  even  though  the  bow-tie 
ideally  is  not  raised  above  the  ground  plane. 

The  Numerical  Electromagnetics  Code  (NEC-2) 
was  implemented  as  the  electromagnetic  simulator 
in  the  GA.  The  planar  structures  under  study  were 
therefore  model^  using  wire  elements.  Figure  2 
shows  an  example  of  the  model  developed  for  the 
BT  antenna.  Ihis  model  essentially  consists  of  a 
best  -  fit  square  mesh  to  which  additional  segments 
are  added  ^ong  die  outside  to  exactly  fit  the  shape 
of  the  antenna.  This  was  necessary  in  order  to 
achieve  the  desired  simulation  accuracy, 


especially  near  the  feed.  The  antenna  was  modeled 
over  an  infinite  ground  plane  and  the  excitation 
was  applied  across  a  segment  of  length  hf  placed 
between  the  bottom  of  the  antenna  and  the  ground 
plane.  A  similar  model  was  developed  for  the  RBT 
by  simply  inverting  the  BT  model.  The 
convergence  of  these  models  was  verified  and  a 
mesh  density  was  chosen  that  was  a  compromise 
between  accuracy  and  complexity  (run-time). 


For  this  study,  only  the  input  impedance  matching 
of  the  antenna  was  optimized  by  Ac  GA,  NEC  was 
used  to  calculate  Ae  input  impedance  at  individual 
frequencies  in  Ac  band  under  optimization.  These 
results  were  used  by  Ae  GA  to  calculate  Ae 
corresponding  VSWR  values  at  each  frequency.  A 
cost  was  Aen  assigned  to  each  design  by  finding 
Ae  maximum  VSWR  value  in  Ae  frequency  range. 
The  GA  sought  to  minimize  Ais  value.  The  GA 
ran  until  Ae  condition  MAX  (VSWR^)  ^  2  was 
met  or  it  was  detnmmed  that  Ae  GA  could  not 
converge  to  a  solution. 

A  steady-state  GA  wiA  binary  tournament 
selection  was  developed.  A  replacement  rate  of 
50%  and  mutation  rates  between  2  -  4  %  were  used 
throughout  The  population  size  used  was  60. 

ra.  RESULTS  /  DISCUSSION 

To  evaluate  Ae  broadband  performance  of  each 
antenna  type,  Ae  G  A  was  run  a  number  times  wi  A 
difTerent  bandwidA  requirements.  For  each  run, 
Ae  bandwidA  was  centered  about  1  GHz.  For 
example,  if  a  fractional  bandwidA  of  20%  was 
desir^  Ae  GA  would  optimize  from  900  MHz  to 
1 100  MHz.  It  was  foimd  that  a  fixed  mesh  density 
of  0.75  cm  could  be  used  to  accurately  model 
antennas  wiA  widely  varying  heights  and  fiare 
angles  in  Ais  frequency  range.  It  was  also  found 
that  these  antennas  could  be  sufficiently 
characterized  by  taking  100  MHz  steps  between 
simulations.  For  Ae  given  example,  Ais  would 
imply  Aat  Ae  OA  would  run  NEC  three  times  in 
or^r  to  evaluate  a  design.  The  overall  dimensions 
of  Ae  antenna  elements  were  nominally  limited  to 
30  cm  X  30  cm  to  keep  GA  run  times  low  and 
because  smaller  antenna  designs  were  desired.  The 
maximum  parameter  values  were  also  limited  to 
H  -  30cm,  a  =  120®,  hf  *=  O.Scm.  A  50fl  feed  line 
was  assumed  Arougfrout  all  simulations.  Figures  3 
and  4  give  Ae  results  of  GA  runs  wiA  increasing 
fractional  bandwidA  requirements  for  Ae  BT  and 
RBT  antennas,  respectively. 

For  the  given  restrictions,  Ae  BT  GA  was  unable 
to  find  a  solution  for  fractional  bandwidths  greater 
than  20%,  while  Ae  RBT  OA  found  solutions  for 
bandwidths  of  up  to  80%.  It  is  clear  from  Aese 
results  Aat  Ae  RBT  is  able  to  achieve  a 
significantly  larger  bandwidA  than  Ae  BT  wiA  a 
dramatically  smaller  size.  Foa-  instance,  Ac  BT 
design  generated  by  Ae  40%  bandwidA  GA  run 
requires  over  5.6  times  as  much  surfece  area  as  Ae 
RBT. 


[ractionai 
bandwidth  (%) 

H 

(cm) 

(deg) 

hf 

(cm) 

cost 

20 

13.33 

95.3 

0.461 

1.90 

30 

14.25 

91.0 

0.416 

2.34 

40 

14.25 

92.2 

0.439 

2.60 

FIGURE  3.  Results  of  GA  optimization  on  the 
bow-tie  antenna. 
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bandwidth  (%) 

H 

(cm). 

a 

(deg,) 

hf 

(cm) 

cost 

20 

9.67 

43.1 

0.37 

1.18 

30 

10.00 

42.4 

0.46 

1.25 

40 

10.33 

38.8 

0.46 

1.34 

50 

10.50 

36.5 

0.46  1 

1.44 

60 

11.00 

35.3 

0.48  1 

1.57 

70 

11.00 

33.0 

0.66  I 

1.71 

80 

11.33 

26.1 

0.78  ! 

1.99 

FIGURE  4.  Results  of  GA  optimization  on  the 
reverse  bow-tie  antexuia. 


As  predicted  by  experimental  data  [5],  even  for 
lower  bandwidA  requirements,  Ae  BT  requires  a 
wide  flare  angle  to  match  to  Ae  relatively  low 
impedance  of  Ae  50Q  feed  line.  Just  to  achieve  a 
20%  bandwidA  requires  nearly  Ae  maximum 
widA  allowed  in  Ae  GA.  Since  no  more  widA  is 
available  to  try  wider  flare  angles  or  taller  heights, 
Ae  BT  is  unable  to  achieve  wider  bandwidAs. 

The  RBT  is  much  more  flexible  in  terms  of  Ae 
matching  quality  vs.  bandwidA  design  trad©K>ff. 
This  is  made  more  clear  from  an  examination  of 
Ae  impedance  plots,  which  are  presented  in  Figure 
5  for  selected  reverse  bow-tie  designs.  For  smaller 
bandwidA  requirements,  Ae  GA  selects  wider  flare 
angles,  which  leads  to  smaller  resistance  values 
wiAin  Ae  optimization  bandwidA.  As  can  be  seen 
in  Figure  6  which  shows  Ae  corresponding  VSWR 
values,  Ae  antenna  is  very  well  matched  at  Ae 
center  of  Ae  band.  To  achieve  a  wider  bandwidA, 
Ae  GA  must  select  a  narrower  flare  angle  to  raise 
Ae  resistance.  This  results  in  a  worse  match  at  Ae 
center  of  Ae  band,  but  a  better  match  overall  so 
Aat  a  wider  bandwidA  is  achieved.  It  is  mteresting 
to  note  that  Ae  height  only  increases  slightly  for 
higher  bandwidA  requirements.  This  is  despite  Ae 
fact  Aat  Ae  current  paA  lengA  has  been  shortened 
due  to  Ae  narrower  flare  angle  and  Ae  lower  band- 
limit  has  been  reduced.  This  implies  Aat  Ae  RBT 
uses  its  frill  surface  area  efficiently  to  achieve 
broadband  performance.  It  is  also  interesting  to 
note  Aat  while  Ae  BT  G  A  sought  roughly  Ae  same 
feed  height  regardless  of  Ae  bandwidA 
requirement,  Ae  vahie  of  Ais  parameter  for  Ae 


RBT  steadily  increased  with  increasing  bandwidth. 
For  the  RBT,  the  feed  height  serves  as  a  means  by 
which  to  optimize  the  matching  of  the  element  to 
the  feed  line.  It  mostly  affects  the  reactance  of  the 
antenna,  scaling  the  reactance  downwards  (more 
capacitive)  for  shorter  feed  heights,  and  upwards 
(more  inactive)  for  taller  feed  heights.  For 
narrower  flare  angles,  the  smaller  element  width 
leads  to  a  reduced  inductance.  Therefore  the  feed 
height  must  be  raised  in  order  to  compensate  for 
the  reduced  elment  inductance  and  achieve 
optimal  impedance  matching. 

The  RBT  designs  generated  by  the  30%  and  70% 
bandwidth  GAs  run  were  both  constructed  out  of 
thin  sheet  copper  and  measured  over  a  ground 
plane  of  dimensions  48  cm  x  48  cm.  These  results 
are  also  shown  in  Figure  6.  The  agreement 
between  the  measured  and  simulated  response  is 
good,  especially  across  the  optimization 
bandwidth. 


FIGURE  5,  Simulated  input  irapedsnee  for  the  30%  end  70% 
bandwidth  reverie  bow>tie  designs. 


IV.  CONCLUSION 

It  has  been  demonstrated  that  the  GA  is  an 
effective  means  of  evaluating  the  broadband 
performance  of  PFMM  antennas.  It  was  shown 
that  the  RBT  antenna  can  achieve  a  much  wider 
impedance  matching  bandwidth  than  the  BT  with 
significantly  reduced  sizes.  The  results  from  die 
GA  illustrate  a  design  metfiodology  for  RBT 
antennas  which  allows  a  tradeoff  between 
matching  quality  and  bandwidth.  This  tradeoff  can 
be  achieved  without  significantly  changing  the 
overall  size  of  the  antenna.  Experimental  results 
confirmed  the  findings  of  the  GA. 
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